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ABSTRACT 


A  computer  model  for  axially  symmetric  dielectric  radomes  in  the  near  field  of 
a  circular  aperture  has  been  developed.  This  model,  based  on  a  method  of 
moments  solution  for  bodies  of  revolution,  is  used  to  determine  the  electric  field 
pattern  of  a  linearly  polarized  circular  aperture  radiating  through  a  dielectric 
radome.  The  program  is  written  in  FORTRAN  and  can  accommodate  rotationally 
symmetric  radome  shapes.  Arbitrary  illumination  distributions  for  the  aperture  can 
be  specified,  and  phased  scanning  of  the  aperture  in  both  azimuth  and  elevation 
is  allowed.  Computational  run  time  has  been  reduced  by  taking  advantage  of 
mode  symmetry.  A  separate  program  has  been  developed  to  compute  gain,  and 
together  these  programs  are  used  to  determine  the  pattern  degradation  and  gain 
loss  due  to  the  presence  of  a  radome  in  the  near  field. 
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I.    INTRODUCTION 

A  radome  is  a  structure  designed  to  protect  an  antenna  in  its  operating 
environment.  Aircraft  radomes  must  protect  radar  antennas  from  high 
aerodynamic  forces  while  minimizing  drag,  weight,  and  interference  on  its 
electrical  operation.  The  presence  of  a  radome  can  affect  the  gain,  beamwidth, 
sidelobe  level,  and  direction  of  boresight,  as  well  as  change  the  VSWR  and 
antenna  noise  temperature.  This  thesis  involves  the  continued  development  of  a 
computer  code  that  predicts  antenna  radiation  patterns  looking  through 
aerodynamic  radomes. 

The  ability  to  accurately  predict  the  effect  of  a  radome  on  the  operation  of  an 
antenna  early  in  the  design  process  is  desirable.  Typically,  designers  of  the 
earliest  radomes  were  concerned  more  with  structural  considerations  than 
electrical  performance.  However,  after  the  near  failure  of  the  U.S.  Army's 
"Pathfinder"  program  of  radar-equipped  B-24  "Liberator"  bombers  due  to  poor 
electrical  performance  of  their  radomes,  it  became  clear  that,"  There  was  a  serious 
need  for  developing  sound  theoretical  procedure  for  the  electrical  design  of 
radomes."  [Ref.  1]  For  most  practical  radomes  the  surface  geometry  is  too 
complex  to  yield  closed-form  mathematical  expressions  describing  their 
performance.  It  is  possible  to  predict  radome  performance  using  numerical 
solutions  of  Maxwell's  equations.  The  program  described  in  this  thesis  is  based 
on  a  numerical  technique  called  the  method  of  moments  (MM).  The  particular 
form  of  the  solution  used  here  is  based  on  a  formulation  by  Mautz  and  Harrington 
[Ref.  7]  for  bodies  of  revolution.  Francis  [Ref.  2]  developed  a  computer  code 
LDBORMM.F  that  specialized  the  Mautz  and  Harrington  solution  to  ogive- 


shaped  radomes.    There  are  no  restrictions  on  the  location  of  the  radome  with 
respect  to  the  antenna,  and  therefore  near-field  radomes  can  be  analyzed. 

The  goal  of  this  thesis  was  to  extend  the  capabilities  of  this  program. 
Significant  improvements  include:  taking  advantage  of  symmetry  to  reduce 
computation  time,  applying  realistic  radar  antenna  patterns  for  comparison  to 
experimental  results,  and  the  computation  of  gain  loss  due  to  the  radome. 

A.    DESCRIPTION  OF  THE  PROBLEM 

The  basic  system  analyzed  is  a  radome  consisting  of  a  body  of  revolution 
(BOR)  located  coaxial  with  a  circular  antenna,  as  in  Figure  1.1.  The  BOR 
constraint  is  satisfied  by  many  aerodynamic  radomes  and  greatly  simplifies  the 
MM  solution  ,  which  will  be  described  in  more  detail  in  following  chapters.  The 


Axis  of 
Symmetry 


Figure  1.1    Radome  and  Antenna 


ogive  shape  of  the  radome  in  Figure  1.1  is  the  primary  one  of  interest,  but  the 
program  can  also  model  other  bodies  of  revolution  including:  a  cone,  disk, 
hemisphere,  and  parabola,  as  shown  in  Figure  1.2.  These  were  primarily  added  for 
validation  purposes. 


There  are  no  constraints  on  the  size,  shape  and  location  of  these  objects 
provided  they  are  coaxial  with  the  antenna.   However,  because  MM  requires  the 


Figure  1.2  Optional  Radome  Shapes 


solution  of  a  matrix  equation,  computer  run  time  and  memory  will  limit  the  size  of 
the  bodies  that  can  be  analyzed.  A  unique  feature  of  this  program  is  that  it  takes 
into  consideration  the  fact  that  the  radome  may  be  in  the  near  field  of  the  antenna 
as  in  Figure  1.3.  The  antenna  aperture  used  in  the  program  is  circular,  but  the 
program  can  be  easily  altered  to  accommodate  other  shapes.  There  are  no 
constraints  on  the  aperture  size,  scan  angle,  or  polarization. 

The  remainder  of  this  thesis  is  divided  into  three  major  sections.  Chapter  II 
covers  theoretical  development  of  the  problem  and  method  of  solution.  The  MM 
solution  is  described  in  Chapter  II  along  with  a  discussion  of  mode  symmetry  and 
how  it  can  be  used  to  advantage  in  obtaining  the  radome  current.  Chapter  II  also 
describes  how  the  total  field  is  calculated  from  the  currents.  Chapter  III  covers 
the    programs    used   for   generating   patterns    based   on    the    methods    of 


Chapter  II.  The  first  part  describes  the  program  RADOME.F  which  calculates  the 
current  and  electric  field.  The  second  part  describes  the  program  GAIN.F  which 


Figure  1.3  Radome  Located  in  Near  Field  of  Antenna 


calculates  the  gain  of  the  system.    The  final  chapter  discusses  the  results  and 
presents  several  conclusions  based  on  the  calculated  radome  patterns. 


II.  ANALYSIS 

A.    GENERAL  DESCRIPTION  OF  ANALYSIS 

The  objective  of  the  analysis  is  to  determine  the  radiation  pattern  of  the 

antenna  and  its  gain  in  the  presence  of  a  radome.   In  order  to  analyze  the  system, 

the  spherical  coordinates  of  Figure  2.1  are  used.    The  total  electric  field  in  the 
direction  of  unit  vectors  (0,0)  is  required  to  compute  the  gain  in  the  far  field.  The 

total  electric  field  is  the  sum  of  the  electric  field  due  to  the  currents  on  the 
antenna  (radiated  field),  and  the  electric  field  due  to  the  currents  on  the  radome 
(scattered  field).  The  procedure  is  to  assume  a  current  distribution  on  the 
antenna,  then  use  the  MM  to  determine  currents  induced  on  the  radome.  In 
practice,  scattering  from  the  radome  back  toward  the  antenna  will  cause  a  change 
in  the  initial  current  distribution.  This  effect  will  be  small  for  a  well-designed 
radome  and  is  ignored  in  the  analysis.  Scattering  back  toward  the  antenna  will 
primarily  affect  the  VSWR  looking  into  its  terminals. 
1.      Determination  of  Gain 

Gain  is  a  measure  ot  the  ability  of  an  antenna  to  concentrate  radiated 
power  in  a  particular  direction.  The  maximum  value  of  gain  for  a  lossless  antenna 
is  the  directivity  Go,  which  is  the  ratio  of  the  maximum  radiation  intensity  to  the 
average  radiation  intensity  [Ref.  3:p  610]. 


4kU  0,0 

Go  =  —4 Znm,  (2-1) 


where  Uy6,o]    is  the  radiation  intensity  and  dQ.  =  S\T\Od6d0  is  the  differential 
solid  angle.  Throughout  this  paper  a  lossless  antenna  will  be  assumed  and  the 


Figure  2.1   Coordinate  System 

terms  gain  and  directivity  used  interchangeably. 

Radiation  intensity  is  the  time-averaged  power  per  unit  solid  angle,  and 
is  proportional  to  the  magnitude  of  the  electric  field  squared 


U  = 


(2-2) 


2ti 


where 


E 

= 

Ea 

+ 

Es 

(2-3) 


and  Ea   is  the  electric  field  from  the  antenna,  Es  is  the  electric  field  scattered  by 
the  currents  induced  on  the  radome,  and  r|  is  the  impedance  of  free  space. 
2.      Determining  the  Electric  Field 

The  total  electric  field  is  calculated  from  the  current  distributions  on  the 
antenna  and  radome  as  illustrated  in  Figure  2.2.  A  current  distribution  Ja  on  the 
antenna  is  assumed  known.    The  currents  induced  on  the  radome  are  dependent 
on  the  current  distribution  on  the  antenna  as  shown  in  Figure  2.3,  and  this  is 
taken  into  account  in  the  method  of  moments  solution. 

Maxwell's  equations,  are  used  to  determine  the  electric  field  at  a  point  in 
space  due  to  a  current  density  J.  When  both  the  radome  and  antenna  are 
present,  the  total  current  density  is  J  =  Ja  +  Js.  A  solution  of  Maxwell's 
equations  is 


E  =  -io)B-^-V(v.fl),  (2-4) 

G)U£      V  ; 


where 


e"jkr 


R  =  JUfj(r,e,<|>)  —  ds\  (2-5) 

4ttjj  r 


and 


r  = 


R-R 


(2-6; 


In  (2-4)  R  is  the  vector  potential,  R   the  position  vector  to  the  observation  point 
and  R'  the  position  vector  to  a  source  point  on  the  antenna. 


E  =  Ea  +  E 


a   .    cs 


Figure  2.2  Electric  Field  from  Currents  on  Antenna  and  Radome 


Ea(r,e,<l)) 


Figure  2.3  Electric  Field  Impinging  on  Radome 


B.    ANALYTICAL  TECHNIQUE 

Closed  form  mathematical  solutions  of  the  gain  and  radiation  patterns  are  only 
possible  for  a  limited  number  of  applications  that  neglect  the  presence  of  the 
radome.  In  general,  determining  the  radiation  pattern  and  gain  with  the  radome 
requires  a  numerical  solution  of  Maxwell's  equations.  The  radome  geometry 
results  in  an  unknown  current  distribution  so  (2-4)  cannot  be  solved  directly. 
However,  boundary  conditions  can  be  applied  on  the  radome  surface  which  yield 
an  integral  equation  for  the  current.  The  integral  equation  will  be  solved 
numerically  using  the  method  of  moments. 

1.      E-field  Integral  Equation  and  the  Method  of  Moments 

The  method  of  moments  (MM)  is  a  procedure  used  to  solve  an  integral 
equation  obtained  by  enforcing  (2-4)  on  the  radome  surface  [Ref.  5].  Equation  (2- 
4)  is  used  to  determine  the  electric  field  Ea  incident  upon  the  surface  of  the 
radome  from  the  currents  on  the  antenna  Ja.  Initially  the  surface  will  be  assumed 
a  perfect  electric  conductor  (PEC),  in  which  case  the  total  E  tangent  to  the  surface 
of  the  radome  is  equal  to  zero.  Later  a  correction  term  will  be  added  to  the  perfect 
conductor  solution  to  account  for  the  dielectric  properties  of  the  radome. 

The  boundary  condition  on  the  PEC  radome  is 

ES,an  =  -Ea,an.  (2-7) 


As  applied  to  (2-4), 


E,ana  =  Mi \\  ±Gds  +  J-v  v.\\  JsGds 

J  J  CO£  , 


(2-8) 


where 


e-Jkr 


G  =  - ,  (2-9) 

4/rr 


and 


r  =  R  -  R 


(2-6) 


Equation  (2-8)  is  the  E-field  integral  equation  (EFIE).  The  only  unknown  is  Js, 
and  once  determined  it  is  possible  to  calculate  the  scattered  field  by  applying  (2-4) 
to  the  radome. 

To  solve  (2-8),  the  current  Js  is  expanded  into  a  series, 

Js=Il,J,-  (2-10) 

1=1 

The  J,   are  basis  functions,  and  the   I,  are  the  expansion  coefficients  to  be 
determined. 

Mautz  and  Harrington  applied  this  technique  to  the  specific  case  of 
surfaces  S  generated  by  revolving  a  plane  curve  about  the  z-axis  [Ref.5].  The 
plane  curve  is  approximated  by  a  series  of  straight-line  segments  connecting  Np 

points.  The  surface  has  a  circular  cross  section  in  the  azimuth  variable  0,  thus  the 
body  is  represented  by  a  series  of  "faceted  rings"  as  shown  in  Figure  2.4. 


10 


The  surface  coordinate  system  is  cylindrical  lp,0,z),  and  t  is  a  vector 

tangent  to  the  surface  and  orthogonal  to  0  at  every  point.  In  order  for  the  series 
expansion  of  the  current  (2-10)  to  approximate  an  arbitrary  surface  current  density 
Js,  independent  sets  of  functions  are  defined  as  follows: 


Triangle  Functions  in  t 


Pulse  Functions  in  Q 


Figure  2.4  Segmented  Radome  with  Current  Basis  Functions 
Depicted 


-T.t 


int:     J'„,  =t 


in  0\     J* .  =  0 

nj 


P. It 


,jno 


,jno 


/7  =  0,±1,  ...,±~        J  =  1,2,  ...,NP-2 


/i  =  0,±1l...,±oo        j  =  ],2,     ,Np-]      (2-11 


Tjtt)  is  a  triangle  function  extending  over  two  segments,  and  Pj[t)  is  a 

pulse  function  extending  over  a  single  segment.    For  an  explanation  of  these 
functions  see  [Ref.  l:pp  16-18].  Each  value  of  n  is  referred  to  as  an  azimuthal 
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mode.    The  series  expansion  of  (2-10)  in  terms  of  the  chosen  basis  functions 
becomes. 


J.=  I 


\p-2  Np-I 

y  /' X  +  y/#  JV 

7=1  7=1 


(2-12) 


In  order  for  the  series  to  be  an  exact  solution,  the  sum  must  be  extended  over 
infinite  modes  -°°</7<°o.  Of  course  the  series  must  be  truncated,  and  the 
maximum  number  of  modes  for  convergence  nmaH  is  dependent  upon  the  antenna 
scan  angle  and  radome  size. 

When  the  series  expansion  is  applied  to  the  integral  equation  (2-8), 


.  +~     Np-2 

el  =nc  \\ 

n=—  7=1  s, 

_  +-    vp-i 

Ea    =  y  y  / 

•-tan  ^      ^  *„ 

n=--  7=1         J, 


M,j' g-^(v'.j'  )vg 


G)£ 


-J 


"7 


JcquJ*G--±  V'.JV  VG 


tfs' 


(2-13) 


In  order  to  solve  for  the  unknowns  l'nJ  and  1°  ,  both  sides  of  (2-13)  are  multiplied 


"7 


by  testing  functions  Ufmi  and  UP   ,  and  integrated.    Using  Galerkin's  Method 

[Ref.  8],  the  testing  functions  are  chosen  to  be  the  complex  conjugates  of  the  basis 
functions  J''", 

nj 


^^  mi 


Uj°     =  0  _LLJ.  e-  7™ 


(2-14) 
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Multiplying  both  sides  of  (2-13)  by  each  of  the  testing  functions  and 
integrating  results  in  a  set  of  2Np  -3  simultaneous  linear  equations  for  each  value 

of  n.  In  matrix  form,  these  equations  are  written  as 


U 


Z    I 


(2-15) 


where 


U 


is  the  excitation  vector. 


is  an  impedance  matrix  and    I   contains  the 

unknown  current  coefficients.     Submatrices  associated  with  the   t,   and    ° 
components  (2-15)  can  be  identified  as  follows: 


V 

n 

= 

"z". 
z" 

n 

z'° 

n 

yoo 
n 

/• 

n 

/I  =  0>±1,±2 /ln 


(2-16) 


The  unknown  coefficients  are  determined  by  solving  the  matrix  equation 


1 

= 

z 

-1 

u 

(2-17) 


The  MM  solution  for  the  unknown  current  coefficients  l'nj  and  1°.  in  vector  form 
is 


n 

z" 

*-n 

yot 
n 

zh 

n 

yoo 
n 

-lr 

'lf.\ 
U' 

n 

(2-18) 


Once  these  are  determined,  it  is  possible  to  compute  the  electric  field  due  to  the 
current  density  distribution  Js  on  the  surface  of  the  radome  at  any  point  in  space. 
The  elements  of  each  submatnx  of  the  excitation  vector   U   are 
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t 


U'    =  J-  ff/Z/'   «Ea   ds 

$ 


(2-19) 


and  the  elements  of  each  of  the  submatiices  of  the  impedance  matrix 


are 


rll 


5,    L 


s,  5/    L 


*,  L 


j^j'„^-^(v.j'n/)v^ 

COE  v  ; 

i^J^-^(V'.Jif-)vC 
0)6:  v 


coe 


J 


jcovJ*G—±V'.J*WG 


coe 


ds  ds 
dsds 
dsds 
dsds. 


(2-20) 


The  fact  that  the  testing  functions  If(J;J  are  orthogonal  in  q  to  J',,',0  due  to 
the  presence  of  the  exponentials  has  been  used  to  eliminate  elements  of  (2-20) 
with  n*m  [Ref.  6].  This  allows  each  mode  n  to  be  treated  independently, 
thereby  simplifying  the  matrix  inversion  problem.  Thus,  by  constraining  the 
surface  to  be  rotationally  symmetric,  the  azimuthal  dependence  of  the  current  is 
expressed  as  a  sum  over  exponentials  (a  Fourier  series)  and  each  mode  n  is 
independent. 

a.     Impedance  of  Thin-Shell  Dielectric  Radomes 


The  impedance  matrix 


of  the  MM  solution  assumes  that  the 
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radome  surface  is  a  PEC.  where  (2-7)  is  applied  as  a  boundary  condition  of  (2-4). 
For  thin-shell  dielectric  radome s.  a  correction  term  is  added  to  the  impedance 
matrix  that  allows  the  modeling  of  radomes  of  arbitrary  dielectric  constant  [Ref.9]. 
When  the  electromagnetic  wave  from  the  antenna  is  incident  on  a 
dielectric  radome  that  has  an  intrinsic  impedance 


=  W  (2-21 


n 


different  from  that  of  free  space  (ri0),  there  will  be  reflection  and  refraction  as 
shown  in  Figure  2.5  [Ref.  3:p  397].  The  total  fields  are 

H    =H$+Ha  (2-22) 

These  fields  satisfy  Maxwell's  equations 

Vx(Ea+EsU-j^o(Ha  +  H$]  (2-23) 


and, 

Vx(Ha  +  HsU  jcoe0[V+V)  +  jo)(e-e0)[V+V),  (2-24) 


where  the  second  term  exists  only  when  there  is  material  present  with  a 
permittivity  different  than  free  space.  This  term  is  a  polarization  current  J.  where 

J  =  jco(e-e0)l.  (2-25) 
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Maxwell's  equations  are  also  satisfied  by  the  field  from  the  antenna    Ea,Ha). 
Since  IE, HI  and  (Ea,Ha)  satisfy  Maxwell's  equations,  by  superposition  and  (2-22) 
Es ,  Hs    also  satisfy  Maxwell's  equations.  Therefore 


VxHs  =  jcoU  +J 


(2-26) 


Figure  2.5  Discontinuity  of  Permittivity  for  Dielectric 
Radomes 


For  a  radome  consisting  of  a  thin  dielectric  shell,  the  polarization 
current  has  a  large  tangential  component  and  small  normal  component.  If  only  the 
tangential  component  is  significant,  then  (2-8)  can  be  modified  to  include 
dielectric  properties  of  the  radome  as  follows: 
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E     "  - 

Man      _ 


jco 


^-  +  jwtipsBds 


coe 


V  .  j]  JsGds 


(2-27) 


The  first  term  on  the  right  side  is  the  impedance  loading  term;  the  integro- 
differential  portion  remains  unchanged  from  the  PEC  assumption. 

After  multiplying  both  sides  of  (2-27)  by  the  testing  functions 


tU^^LUmi)  and  integrating,  the  impedance  matrix 


of  (2-15)  becomes 


|Z]  =  [ZMM]  +  [ZL], 


(2-28) 


where  the  elements  of 


-MM 


are  given  by  (2-20),  and 


Z" 

■jto 

*-ln 

■jot 
*-Ln 


J'J 


J'V 


J'J 


=  \\U!!.JsjnZLds 

s 

=  [[  UJf  •  J°sjnZLdS  =  0 ,  (by  orthogonality  of  t,  and  o  ) 

s 

=  jj  Iff  •  J'sjnZLds  =  0,  (by  orthogonality  of  0,  and  t) 


l?a]u=\\W!'*sjnzLds, 


(2-29) 


and 


Zl  = 


JQ)[£-£0)t 


(2-30) 


The  thickness  t  in  the  impedance  loading  term  is  necessary  to  convert  the  volume 
current  density  J,  of  Figure  2.5  to  a  surface  current  density  Js. 
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Solving  (2-17)  yields 


Zmm+Z, 


u 


(2-31) 


The  radome  impedance  ZL  can  be  rewritten  in  terms  of  dielectric  constant  £r, 
thickness  to  wavelength  ratio  n, .  and  the  loss  tangent  tan  8  [Ref.  2:pp  45-46], 


Zl 


60 


n 


j[er  -1)  +  ef  tan 8 


(2-32) 


Note  that  both  (2-29)  and  (2-32)  allow   the  modeling  of  radomes  whose 
permittivity  may  be  complex,  e  =  e'  -  je"  . 
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2.       Mode  Symmetry 

When  the  surface  current  density  Js  is  expanded  into  a  series  as  in  (2- 
12),  the  sum  is  extended  over  aznnuthal  modes  n  =  0,±1,±2, ■•■,±f?ma„.    The 

mode  index  appears  in  the  e\ponential  term  of  basis  functions  J'n°    (2-11)  and 

testing  functions  LU'°   (2-14).  The  exponential  factor  can  be  written  as, 

ejn  =  cosine)  +  jsinino).  (2-33) 

For  each  n  there  is  a  current  \  ector  l„  given  by  (2-18).  Mode  symmetry  refers  to 
a  relationship  between  the  positive  and  negative  mode  vectors  l„  and  l_n.  If  mode 
symmetry  exists,  then  only  one  matrix  equation  needs  to  be  evaluated  rather  than 
two.  This  greatly  reduces  the  computational  effort. 

Recall  the  definition  of  the  excitation  vector   U 
If.,  =-\\Ufm.vJds 

s 

The  electric  field  from  the  antenna  is 

Ea  =  EajR,6,o)fl  + E'JR, d, 0)6 +  Ea,(R, 6,0)0,  (2-34) 


which  on  the  radome  surface  can  be  expressed  in  terms  of  the  two  orthogonal 
components  t  and  0 

V  =E?(t,0)i  +  EaJt,<l>)0.  (2-35) 
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In  the  far  field  of  the  antenna,  the  R  component  in  (2-34)  decays  to  zero  and  Ef 
consists  only  of  E*.  Therefore  as  shown  in  Figure  2.6, 


E?(t,o)  =  E?(t)cos0 
EaAt,o)  =  Ea(t)sin0. 


(2-36) 


Cut  through  sphere  in  0=  constant  plane 


Radome  Curve 


Figure  2.6  Components  of  the  Antenna  Electric-Field  on  the  Radome 


Plugging  these  far-field  components  into  (2-19) 


U'm  =  J  TJ^)p(t)p(t)dt  J  cos0[cos(n0)  +  jsin(n0)]t/0 


f=8 

I 


o=8 
2k 


U">  =  1  Pj^P[t)p[t)ilt  1  sin0[cos(n0)  +  jsin(/?0)]r/0. 

f=8  Ps^   ' 


o  =  8 


(2-37) 
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Only  the   o  integration  is  dependent  on  n     Expanding  the  integral  in  o  and 
explicitly  including  the  sign  ol  n  gives 


J  coso\cos[±no]j  +  jsin{±no)\j0  = 

=  9 

lit  2-T 

J  cosocos[±n0)d(D  +  j  Jcos0sin(±no)(/o,  (2-38) 


o  =  B 
2/r 


o  =  9 


where 


2- 


j  cos<z)Sin[±/?<2>]£/0  =  0 


(2-39) 


o  =  0 


due  to  the  fact  that  the  integrand  is  an  odd  function  of  0,  and  the  integration  is 
over  one  period  0  <  (f>  <  In.   By  trigonometric  identity 


2-t  2rz  2/r 

J  cos0cos(±f?0)£/tf>  =  j  cos  i±n-]\0m+  J  cos  (±n  +  \U 

0=9  =9  «=9 


do     (2-40) 


The  right  hand  side  is  only  non  zero  if  it  =  ±1.  Furthermore,  because  cosine  is  an 
even  function,  the  term  U*ni  is  independent  of  the  sign  of  n.  Thus 


u(nj-^+nj 


(2-41) 


For  the  second  equation  of  (2-37 ).  the  sign  of  U°   is  dependent  on  the  integral 
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where 


2k 


\  sin 


o 


cos(±A70  +  jsin  [±n0 


2k 


]t/0  = 


2*  2-t 

J  smocos[±no]jdo  +  j  Jsin0sin(±/?0)(/0,  (2-42) 


0  =  8 


C  ■ 

jsinocos(±nc))(/0  =  0, 


o  =  0 


by  symmetry.  By  tiigonometnc  identity 


12-43) 


2,t  2t  2^- 

J  sin0sin[±/70)d0  =   )'  cos  (tn-IJoW-  f  cos 

o=9  o=B  o=B 


±n  +  \)0 


dp      (2-44) 


As  can  be  seen  in  (2-44),  the  excitation  vector  U'   is  dependent  on  the  sign  of  n 
due  to  the  presence  of  the  negative  sign.  Therefore 


U°     =  -U° 

-nj  +nj  • 


(2-45) 


In  the  near-field  of  the  antenna  the  0  dependence  is  more  complicated 

than  the  sine  and  cosine  functions  in  (2-36).    However,  it  has  been  verified 
numerically  that  E*  and  E J  are  always  even  functions  of  0  and  E '  is  always  odd. 

This  conclusion  assumes  that  a  symmetric  amplitude  and  linear  phase  distribution 

excites  the  antenna. 

In  light  of  the  above  results,  the  elements  of  the  impedance  matrix   ZMM 

of  (2-18)  have  the  following  relationship  between  positive  and  negative  modes: 

[Ref.  6] 
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-MM 


[ZMM,-n;j;j        |^M(-fl)jy 


yo  t 
*-MM(-n) 


J 'J 


J<r>Q 

tMM(-n) 


'7 


■MM(  +  ff 

tMM(+»; 


U 


-Z 


to 

MMlt/l 


■joo 


'J 


'J    _ 


(2-46) 


By  inspection  of  (2-30),  the  elements  of  the  load  impedance  matrix  ZL  are 
independent  of  the  mode  /7.  because  W\°n  and  J'j^  are  complex  conjugates  of  each 
other.  Therefore, 


zL 

L  L 

r 

)_ 

U 

0 

z" 

)_ 

»J 

0 

0 

■j  00 

_ 

U  . 

0 

■j  00 
.    L(+" 

_ 

t/_ 

(2-47) 


Rewriting  (2-31)  results  in 


U 


2-48! 


where 


ZMM   +  ZL 


n-1 


(2-49) 


By  taking  the  inverse  through  partitioning  [Ref.  8]  it  can  be  shown  that  the  mode 
symmetry  survives  the  inversion,  and 


Y" 

■  -nij 

u   ' 
_L    ~n,J  . 

uh 

-nij 

y; 

-nij 

.' 

-i-n/j 

_uof 
+  nij 

_uio 

+  nij 

u 

+  nij 

(2-50) 


By  multiplying  out  (2-48), 
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I'  =  V"  u'  +  vh  u° 

*n         ■  mjunj  T  *  nij     nj 


I    =  ¥'.!/'    +VM.[/V.  (2-51) 


Equations  (2-41 ).  (2-45),  and  (2-50)  relate  the  positive  and  negative  modes  of  (2- 
51),  as  follows: 


I'    =v"nilu'ni  -¥"  ..(-[/•  .  =+i' 

-n  -nij     *nj  +n/j  \         +nj  j  +n 

V    =  -¥'    U[n , +  ¥T \-U<  .\=-l*  .  (2-52) 

-n  -nij     +nj  +mj  \         +nj  +n  y  ' 


This  is  an  extremely  important  result.   Using  mode  symmetry  reduces  the  number 
of  calculations  required  to  find  the  current  coefficients  I''*  by  nearly  half  in  most 

cases. 
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3.      Measurement  Matrices  for  the  Scattered  Field 

Once  the  current  Js  is  obtained  the  scattered  field  Es  can  be  determined 
by  summing  the  contributions  of  the  individual  current  elements  in  the  far-field  as 
shown  in  Figure  2.7.  This  is  accomplished  by  use  of  a  measurement  vector  with 
elements 


K  =  \\[n  +  e  +  o)e-Jkr  *dsds. 


(2-53 


In 

far- fie Id 


Figure  2.7  Scattered  Field 


In  (2-53),  a  unit  radiated  plane  wave  is  weighted  by  the  current  element  at  each 
point  on  the  surface.  Neglecting  the  R  component  in  the  far-field,  the 
measurement  vector  separated  into  0  and  0  components  is 


R9m=.lJe-Jk(Ka+yu+a"  ]\ 
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»:-JJe--w"*+-,(^.31)rf«'.  (2-54) 

$ 

where  U,  U ,  and  LU  are  the  direction  cosines, 

u  =  sin#cosc5 
u  =  sin#sin<z> 
iv  =  cosd  (2-55) 

Thus,  the  components  of  the  scattered  field  at  a  far-field  point  are  the 
superposition  of  all  surface  segments 


4.      Closed-Form  Far-Field  Antenna  Pattern 

When  computing  the  antenna  contribution  to  the  total  field  in  the  far 
zone  the  higher  order  j  terms  can  be  neglected,  and  the  radiation  integral  can  be 

reduced  to  a  closed  form.   The  source  of  this  field  is  assumed  to  be  a  uniformly- 
excited  circular-aperture  antenna  scanned  to  a  direction   [9S,0S).    Uniformly 


excited  refers  to  a  constant  amplitude  and  a  linear  phase  distribution  of  current  on 
the  antenna.  The  /?,  6  and  0  components  in  the  far-field  are 


f;=8 
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-  ilf  p'Jkr  i  s 

where  L0,  LH    Ny   and  N0    are  given  by 

L  =  \\Me-Jkr  C0SU  rfs'  =  0,  (M  =  8), 

N  =  JJ  Je"JJfr  C0SV  f/s'  =  \\[kJh  +\}Jy  +zJzyjkrcosv  ds'         (2-58) 


s 


where  r'  and  y/'  are  as  shown  in  Figure  2.8  [Ref.  10:pp.  455].    If  the  antenna  is 
x-polarized,  then  Ny  -  Nz  =  0  and 

NH  =  \\jHe~Jkrcosv'ds'.  (2-59) 


The  equivalent  in  spherical  coordinates  is 

Ne  =  \\JH  cos  e cos  0e-jkr'zosv'ds' 
N0  =\\jHs\n0eJkrcosv'ds'.  (2-60) 


27 


The  phase  of  the  source  point  relative  to  the  field  point,  r'  COS  y/\  is  a  length  that 
can  be  written  in  Cartesian  coordinates  as 


r' cos  i//  =  h'u  +  y'u  +  z'w 


(2-61) 


Source  Point 


Figure  2.8  Circular  Aperture  Antenna  Source  Point 


where  U,  U,  and  Hi  are  the  direction  cosines  (2-55).   Including  beam  scanning 
adds  an  exponential  phase  factor  to  the  current 


JM  =  J0e 


jk\KV,+yv,+zwt 


(2-62) 
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where  Us.  us  and  Ws  are  the  direction  cosines  of  the  scan  angles  ( 6S,  os\  which 

also  obey  (2-55).    J0  is  the  magnitude  distribution,  assumed  to  be  constant,  and 
[#',*/',  Z' )  are  the  coordinates  of  a  source  point.   For  an  antenna  lying  in  the  x-y 

plane,  Z'  =  0.  Therefore,  NH   can  be  rewritten  as 


Nti  =  J0  cos  h  cos  offe 

s 

=  J0cost>coso[[e 


jk\ h  us  +y  u, )    -jk[HV+y  V )., 


jk{n-lus-u)+y(u,-u] 


ds' 


(2-63) 


$ 


By  setting  H'  -  p'CQS0'  and  y'  =  p'Sind',  the  exponent  term  of  (2-63)  become 


kP' 


cos o[us  -u)  +  s\r\0' [us  -u 


(2-64) 


Defining  the  variables 


R  =  [uc  -u 


B  =  (u,-ul 


(2-65) 


and  considering  these  to  be  orthogonal  vectors  as  in  Figure  2.9,  it  is  possible  to 
rewrite  (2-64)  as 


kp'^R2  +  B2 


coso  +  sin  0-t- 


Vfl2  +  B2 


V/r  +  B2  \ 


2-66! 
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From  Figure  2.9  the  factors  containing  R  and  B  in  the  brackets  can  be  replaced 
with  COS/  and  Sin/ 


kp'4n2  +  fl2[coso  cos /  + sin 0' sin/ 


(2-67) 


By  trigonometric  identity. 


k,y^R2  +  B2 


cos w  -y 


2-68) 


B 


R 


Figure  2.9  Representation  of  R  and  B  as  vectors 


Using  this  form  of  the  exponent  in  the  integral  of  (2-63),  gives 


2,t      a 


Nd  =  J0cos#cosc>  J     \e 


*o-W+''"4'-r)p,dp,d0, 


(2-69) 


o'  =  B  p'  =  8 
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which  conforms  to  the  definition  of  a  Bessel  function.  The  same  integration  is 
performed  in  the  second  equation  of  (2-60),  and  the  same  Bessel  function  results. 
Finally, 


%  =  J0  COS  tiCOSO  ntf 


JAka\R2  +  B2 


2        n2 


ka^R2  +  B 


N0  =  -J0  sino  nal 


2    ,    n2 


J,  kaiR2  +  B 


2    ,    D2 


ka\R2  +B 


2-70; 


where  J,  is  a  Bessel  function  of  the  first  order.    The  far-field  equations  (2-57 
become 


ikn  &~^r 
Ee=-jo±7-L-      -cos#cos<z>  na2 


4/r     r 


jkr 


n2    ,     n2 


2  J,  ka^R2  +  B 


2         n2 


ka^R2  +  B 


r       ,  jkn  e 

4/r     r 


2J,  JraV/r  +  B 


|2    ,     n2 


ka^R2  +  B2 


(2-71) 


where 


r  =  us  -  u  =  sin  ds  cos  os  -  sin  d  cos  <p 
B  =  us  -u  =  s\nds  sin0s  -sin0sin<z>, 


2-72) 


and 
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Rl  +  Bl  =  Jsirr  #  +  sirr  6S  -2sin#$  sin d cost 0-0, J.        (2-73) 


Using  equations  (2-71)   through   (2-73)   are  computationally   less 
intensive  than  a  numerical  integration.   These  equations  were  employed  for  both 
the  far  field  pattern  and  gain  computations. 
5.      Near  Field  of  the  Antenna 

In  the  previous  section,  the  far- field  of  a  circular  aperture  antenna  was 
determined.  When  a  radome  is  located  in  the  near-field,  the  integrals  cannot  be 
reduced  to  closed  form  and  they  must  be  evaluated  numerically.  The  near  field  is 
generally  considered  to  be  the  region  surrounding  the  antenna  where  r  «  \ ,  or 
kr  «  1  [Ref.  10:pp.  106]. 

The  complete  expressions  for  the  near-field  of  a  circular  aperture 
antenna  excited  by  a  current  polarized  in  the  x-direction  are 


E-=^t\\  G,Jh+(h-x')2G2Jh  remedy 

f«=^JJ[(»-"')(y-i/')eA]e-J,;'rf«-dy' 

s' 

Ez=^t\\[[x-H')(z-z')G2JHy^dH'dy\ 


(2-74) 


where 


Gi 


klrl  - 1  -  jkr 


2-2 


G2  = 


3  +  3  jkr  -  k'r 


(2-75) 
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and 


r  =  ^(x-K')  +\y-y)  +[z-z')  .  (2-76) 

The  primed  coordinates  designate  source  points  and  the  unprimed  an  observation 

point  on  the  radome  as  in  Figure  2.10.    [Ref.  11]    S'  is  the  area  of  the  antenna 
aperture.    Since  the  radome  is  rotationally  symmetric,  the  components  lEHEyEz) 

can  be  described  in  terms  of  cylindrical  coordinates  using  the  transformations 


P [H,y, Z   or  PiR,  6,0 


Figure  2.10  Rectangular  Coordinates  for  Circular  Aperture  Antenna 


X'  =  p'COS0' 

W  =  p  sin0' 

dH'dy'  =  p'dp'dQ'. 


(2-77) 


For  an  antenna  scanned  to  [ds,os\  the  current  distribution  is 
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J,  =  «/,p'fl 


jk\x'cos \o'-0   | sin ^  +  f/'sin|  0' -  0    sine 


.    (2-78) 


Applying  (2-75),  2-76)  and  (2-77)  to  equations  (2-74)  results  in 


2/r      a 


f»=<£  J    )J0(p')e*'CO5l,-)sl"e,e-*x 


0  =a  p  =9 


G,  +  (pcos0-p'cos0')  G 


p'dp'do' 


2k      a 


J    K(p' 


jkp'cos[o-9,)s\na,      jig. 


s rJ"  x 


0=8  p  =e 


pcos<z>-p'cos0'  psin^-p'sinp'  G2  p'dp'do 


2k      a 


£z    =^L    J       Jjo(p'jeJ*P'«:os(«-«f|sinefe-j»:r  x 


0  =8  p  =8 


p cos 0-p' cos <f>')zG2  p'dp'do' 


(2-79) 


where 


»'a#-/,1  +  z"z'    +  4pp'sin 


2  0- 


(2-80) 


Since  the  expression  for  the  excitation  vector  requires  spherical  components, 
lEHEyE2  \  are  converted  using  a  coordinate  transformation 


34 


ER  =  EH  s\nn  cos  o  +  E^  sine  sin  o  +  Ez  cos  6 
fw  =  EH  cosh  cos  o  +  Ey  cose  s\no-Ez  sine 

E0  =  -EH  sin  $  +  Ey  cos  o.  (2-81) 

These  are  the  equations  programmed  to  compute  the  E-field  impinging 
on  the  radome.  The  magnitude  of  the  current  distribution  is  allowed  to  vary  as  a 
function  of  p' .  Various  antenna  sidelobe  levels  can  be  modeled  by  simply 
changing  the  function  Ja(p') 
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III.     PROGRAMS 

The  mathematical  formulas  developed  in  the  previous  chapters  were 
computer  programmed  using  FORTRAN.  The  pattern  and  gain  calculations  are 
done  in  two  separate  programs:  RADOME.F  and  GAIN.F,  both  of  which  run  on  a 
Sun  SPARCstation™  under  UNIX.  A  brief  description  of  each  follows. 

A.    DESCRIPTION  OF  THE  PROGRAM  RADOME.F 

This  program,  developed  by  Francis  [Ref.  2],  determines  the  series  expansion 
coefficients  of  the  current  Js  induced  on  the  surface  of  the  radome.  It  also 
calculates  the  total  electric  field  surrounding  the  system,  and  generates  radiation 
patterns  using  MATLAB™.    Originally  called  LDBOR.F,  (an  abbreviation  of 

"Loaded  Bodies  of  Revolution")  in  Ref.  2,  it  has  been  modified  to  take  advantage 
of  symmetry,  and  also  allows  the  antenna  mainbeam  to  scan  in  both  6   and  o  . 

Figure  3.1  is  a  flow  chart  of  the  program.  To  run  the  program,  data  is  entered 
interactively,  except  for  a  data  file  "gaus/gaus(ITH)M  containing  the  integration 
constants  for  Gauss  Quadrature,  which  is  read  in  the  program.  All  integrals  are 
evaluated  using  this  technique  (see  Appendix  A.).  Table  I  is  a  sample  input. 

After  all  data  is  entered,  the  coordinate  points  of  the  surface  are  computed 

and  stored  in  the  arrays  RH(  =  p )  and  ZH(  =  Z).   The  surface  impedance  of  each 

segment  is  stored  in  ZLO.    The  elements  of  the  MM  impedance  matrix  are 

computed  in  subroutine  ZMAT.  and  the  excitation  vector  elements  in  GENEX. 
These  provide    ZMM    and    U    in  equation  (2-17).    Subroutines  DECOMP  and 

SOLVE  solve  the  matrix  equation  using  Gaussian  elimination.  Subroutine 
PLANE  calculates  the  measurement  matrices,  which  are  subsequently  used  to 
compute  the  scattered  field  from  the  radome. 
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Two  antenna  radiation  subroutines  are  included.  CIRCRTP  calculates  the 
electric  field  components  EH,E^.  and  E0  according  to  equations  (2-81 ).    It  is  used 

by  subroutine  GENEX  to  obtain  the  excitation  vector.  The  second  antenna 
subroutine  ANTFF  assumes  that  the  observation  point  is  in  the  far-field.  This  is 
used  in  the  main  program  to  compute  the  antenna  contribution  to  the  radiation 

pattern  E'  in  (2-3). 

Output  includes  an  ASCII  file  "outldbor"  which  lists  all  the  input  parameters, 
radome  geometry  and  impedance,  and  the  pattern  magnitude.  The  current 
coefficients  are  written  on  a  file  "curcoefsdat"  so  that  further  patterns  can  be 
generated  without  recomputing  the  coefficients.  Finally,  MATLAB™  formatted 
files  are  generated  for  plotting  the  6  and  o  components  of  the  elecuic  field. 
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TABLE  I.  SAMPLE  INPUT  RADOME.F 


ENTER  A  LETTER  TO  INDICATE  WHICH  ITERATION  THIS  IS      A 
SELECT  BOR  GEOMETRY  BY  NUMBER 


vlBER 

BOR 

1 

OGIVE 

2 

SPHERE 

3 

CONE 

4 

DISK 

5 

PARABOLA 

ENTER  SURFACE  CURVATURE  (wavelengths)  7.5 
ENTER  ZPRIME.WHERE  CURVATURE  STARTS  (wavelengths)    0 

ENTER  BASE  RADIUS  (wavelengths)  3 
ENTER  FILENAMES  gaus### 

ENTER  THE  FILENAME  IN  T  (NT)  gaus2 

ENTER  THE  FILENAME  IN  PHI  (NPHI)  gaus48 

ENTER  THE  FILENAME  IN  X  AND  Y  gaus20 

ENTER  THE  PLOTTING  INCREMENT  IN  DEGREES  1 

ENTER  THE  HIGHEST  MODE  20 

ENTER  PHI  (observation)  IN  DEGREES  90 

ENTER  SCAN  ANGLES  IN  DEGREES:  THETA,  PHI  30,90 

ENTER  COMPLEX  IMPEDANCE.  OHMS:  (Real,  Imag)  (0,-1700) 

ENTER  THE  ANTENNA  RADIUS  (wavelengths) 1 
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Output  written  to  files  appended  by  IT 


l=Ogive 


i 


2=Spliere 


3=Cone 


4=Disk 


5=Parabola 


Call  OGIVE 
Determine  NP, 
ZH.RH.BASE 
RS.ZP 


Call  CONE 
Determine  NP, 
ZH.RH.BASE, 
RS.ZP 


Call  SPHERE 
Determine  NP, 
ZH.RH.BASE. 
RS.ZP 


Call  PARAB 
Determine  NP, 
ZH.RH.DM. 
FOD.ZP 


Call  DISK 
Determine  NP. 
ZH.RH.BASE 
RS.ZP 


NP=  #  points  on  curve  generating  BOR 

ZH.RH=Coordinates  of  BOR 

DM=Diameter  of  Parabolic  Radome 

FOD=F/D 

ZP=Location  of  antenna  z' 

B  ASE=Base  diameter  of  antenna 


© 


Enter  filenames  for  Gauss  in  t 

and  4) 


Figure  3.1   RADOME.F  Flow  Diagram 
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0 


Weights  and  Abscissas  are  read  in  from 
files  named  Gaus/Gaus(  ITH) 


A 


Read  XT.XP.X.A.AT.AP 
^  ~ 


7 


Inter  DT,  Plotting  Increment 


No 


A 


Enter  Highest  Mode.  MODES. 


Calculate  NBLOCK.NROW 


NBLOC  partitions  current  vector  [I] 

into  by  mode  n 
NROW=total  #  of  individual  elements 

of  [I] 


Read  ITH.SELECTION.BASE 
NP.MT.MPN.NROWS, 
MODES.  NB  LOCK.THS.  1  'HS . 
ARAD.PHIJMP.ZP.RS.R  H. 
ZH,ZLO  from  "curcoefsdai(ITH)" 


Calculate  1T.MHI.MP.MT.N 


Prepares  for  E-field  calculation 


Figure  3.1  RADOME.F  Flow  Diagram    (Continued) 
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0 


Exit  i * »  Loop 

—([X)S2I=1,NPW 


Calculate  ZHB.RHB.ZLOd) 


Assign  surface  impedance  to  each 
segment  of  radome 


No 


Call  ZLOAD 
returns  ZL 


Exit 


F^- 


Subroutine  ZLOAD  calculates 
elements  of  [ZjJ 


DC)4(K)M=1.MHI 


■vLoop 


n=+M 


T^-\ 


Positive  mode  n 


Call  ZMAT 
returns  Z 


IZH 


Subroutine  ZMAT  calculates 
elements  of  [Zj^^j] 


Call  ZTOT 
returns  Z 


Subroutine  ZTOT  calculates 
elements  of  [Z]=[Z  MM  Wzl1 


Call  GENEX 
returns  B 


s 


Subroutine  GENEX  calculates 
elements  of  excitaion  vector  [V" 


Call  DECOMP 
returns  Z 


6 


Subroutine  DECOMP  calculates 
■1 


elements  of  [ZjJ" 


6  : 

Figure  3.1  RADOME. F  Flow  Diagram  (Continued) 
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o© 


9 

A 


Call  SOLVE 
returns  C 


Subroutine  SOLVE  calculates 
elements  of  [I ]=[Z]'1[V] 


Store  CT,  (+n)  current 
elements 


Negative  mode  -n 


Call  GENEX 
Call  ZMAT 
Call  ZTOT 
Call  DECOMP 
Call  SOLVE 


Store  CT,  (-n)  current 
elements 


Write  I  TH. SELECTION  BASE.NP.MT.MP.N 
NROW.MODES.NBLOCK.THS.PHS.ARAD, 
PHI.IMP.ZP.RS.IFF.RH(L).ZH(L),ZL(0),ZLO(L 
GNT.C.NPHI.CGP  to  curcoefsdat(ITH) 


These  parameters  are  to  be  used 
in  Gain  computation 


X     1 

O0 


Figure  3.1   RADOME.F  Flow  Diagram   (Continued) 
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0 


Exit 


DO  500 


=i.it) 


Loop 


Compute  E-field  tor  each  angle 
of  plotting  increment 


Initialize  ET1.ET2 
ppi  pp"> 


Exit 


R 


DO  30()M  =  1.MHI 


Loop 


© 


Compute  Scattered  Field,  Es  (2-13) 


Calculate  EXP1.EXP2 


e+/-jn<j>(2-56) 


Call  PLANE 
returns  R 


Exit 


Subroutine  PLANE  calculates 
elements  of  Measurement  Vector 
R  (2-54) 


DO402L=l,MT 


\Loop 


Calculate  ET1.EP1 


ET1.EP1  are  elements 
of  ES.  (2-56) 


Exit, 


DO  403  L= 


.Loop 

L=MT+1,N   ) 1 


Calculate  ET2,EP2 


ET2.EP2  are  elements 
ofEs,(2-56) 


*L 

0 

Figure  3.1  RADOME.F  Flow  Diagram   (  Continued) 
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0    0 


Call  ANTFF 
returns  ETF.  EPF 


Subroutine 
ANTF1' 
calculates 
tar-field 
of  antenna 


Call  CIRCRTP  returns 
CIRCR,  CIRC.  CIRCP 


Subroutine  CIRCRTP 
claculates  near  field  of 
antenna 


Calculate  ETF.EPF 

from  CIRCR.  CIRCT.CIRCP 


Calculate  ETHF.  EPHF,  ET,  EP, 
ETSCA  T.  EPSCAT.  EC.  EX, 
ECV.  EXV,  ECR.  ECI,  EXR, 
EXI.  PHC.  PHX.  ANG.  ECX 


(8:<J))  components  of  E-field.  total  field, 
scattered  field,  copolarized,  cross  polarized, 
phase,  angle 


Write  E-field  components  to 
separate  files 


E-field  files  are  formatted  for  input  to 
MATLAB^M  for  plotting  purposes 


(       End       ) 


Figure  3.1   RADOME.F  Flow  Diagram  (Continued) 
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B.    GAIN.F 

This  program  determines  the  gain  of  the  system  by  solving  (2-1 ).  In  order  to 
perform  the  integration,  the  current  coefficients  from  RADOME.F  must  be 
provided  so  that  the  radome  scattered  field  can  be  computed.    Input  consists  of 
specifying  an  iteration  letter,  and  selecting  the  number  of  divisions  in  the  6  and  0 
integration.   Input  data  is  then  read  from  a  file  "curcoefsdat"  that  is  appended 
with  the  iteration  letter  chosen. 

Figure  3.2  is  a  flow  chart  of  the  program,  and  Table  II  is  a  sample  input. 
GAIN.F  contains  several  of  the  same  subroutines  that  occur  in  RADOME.F.  They 
include  PLANE,  CIRCRTP,  and  ANTFF. 

TABLE  [I.  SAMPLE  INPUT  GAIN.F 


ENTER  A  LETTER  TO  INDICATE  WHICH  ITERATION  THIS  IS  A 
ENTER  NUMBER  OF  DIVISIONS  IN  THETA  DIRECTION, 

NDIVT  =  ?  2 

ENTER  NUMBER  OF  DIVISIONS  IN  PHI  DIRECTION,  NDIVP  =?  3 
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(     Begin     ) 


Assiiin  PiixiUTieters 


Read  ITH.SI  LECTION.BASE 
NP.MT.MP.N.NROWS, 
MODES.NBLOfK.THS.PHS, 
ARAD.PHI.IMP.ZP.RS.RH, 
ZH.ZLO  from  "curvoefsdatdTH)" 


data  was  written  to  "curcoefsdat(ITH)" 
from  RADOME.F 


I   ReadCKL)    / 


CT(L)  are  elements  of  current  coefficient 
vector  [I]  written  to  curcoefsdat 


L 


Read  XT(I),  AT(  I).  Xtl).  Ad).  XP(I).  AP(I) 


A 


7 


A(I)  are  weights,  and  X(I)  are  adscissas 
which  are  read  from  2aus/saus# 


Enter  NDIVT.  NDIVP 


i 


7 


NDIVT,  and  NDIVP  are  the  number  of 
divisions  in  9  and  0  respectively.  This 
allows  for  easy  variation  of  the  number 
of  integration  steps  without  recomputing 
weights  and  abscissas  for  Gauss 
integration. 


Calculate  SO),  b  integration  angles 
andQ(I).  0  integration  angles 


I 


Set  SUM,  SUMNONE.  UMAX.  UMAXNONE, 
EMAX,  EMAXN(  »NE  to  0 


Figure  3.2  GAIN.F  Flow  Diagram 
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Exit 


11  y  ■ y    1    IHI   > 

■^  Do3(XljJ=l,NDIVP /J*    — 


Summation  of  divisions  in  0 


AXip 


Do  300  .1=1  .NPM I       Vy    ~~    Summation  within  each  division  in  o 

Loop 


Summation  of  divisions  in  8 


Summation 
within  each 
division  in  £ 


E-field  determined  at 
each  THR  ( 9),  and 
PHR  (0)  angle 


Set  ET,  EP.  ET1,  EP1.  ET2.  EP2  to  0 


Exit 


0. 


{  Do305M=l.MHI 


Loop 


Mode,  n.  summation 
of  E-field  (2-13) 


Calculate  EXP1.EXP2 


e+/-jn<J)(2-56) 


Call  PLANE 
returns  R 


V— "x  /""■'x  Subroutine  PLANE  calculates 

elemenus  of  Measurement  Vector 
R  (2-54) 


0  00 


© 


Figure  3.2  GAIN.F  Flow  Diagram   (Continued) 
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0© 


o 


DO  306  L=  1,1 


Calculate  ET1.EP1 


ET1.EP1  are  elements 
of  Es,  (2-56) 


ET2.EP2  are  elements 
of  Es.  (2-56) 


Call  ANTFF 
returns  ETF.  1-PF 


Subroutine 
ANTFF 
calculates 
far-field 
of  antenna 


Call  CTRCRTP  returns 
CIRCR.  CIRC.  CIRCP 


Subroutine  CIRCRTP 
claculates  near  field  of 
antenna 


Calculate  ETF.EPF 

from  CIRCR.  CIRCT.CIRCP 


Calculate  ETCOMB.  EPCOMB 
EMAG,  EMAGNONE.  L1MAG. 
IIMAGNONE.  SUM.  SUMNONE 


j 


Total  E-field  in  9  and  <J>,  magnitude  of  E-field 
with  and  without  radome,  magnitude  of  field 
intensity  with  and  without  radome.  current 
sum  with  and  without  radome. 


0 


Figure  3.2  GAIN.F  Flow  Diagram  (  Continued) 
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Calculate  PR  AD.  PRADNONE 


PRAD=  radiated  power  with 

radome 
PRADNONE=  radiated  power 
without  radome 


Calculate  GAIN.  GAINN(  )NE. 
GDB.GDBNONE 


GAIN=  gain  with  radome 
GAINNONE=  gain  without 

radome 
GDB=gain  with  radome,  (dB) 
GDBNONE=gain  without 
radome.  (dB) 


Write  gain  values  to  "gainout(ITH).m' 


"gainout(ITH).m"  is  formatted  for 

input  to  MATLAB™ 
for  plotting  purposes 


(       End       ) 


Figure  3.2  GA/N.F  Flow  Diagram  (  Continued) 
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IV.     RESULTS  AND  CONCLUSIONS 

A.    VALIDATION  AND  CALCULATED  RESULTS 

Part  of  the  continued  development  of  this  research  is  validation  of  the 
programs.  Test  cases  were  run  to  see  how  the  programs  performed  against 
known  results.  Test  radomes  ranged  from  nearly  nonexistent  disks  located  at 
large  distances,  to  perfectly  conducting  parabolas  located  in  the  near-field.  E- 
plane  and  H-plane  patterns  defined  in  Figure  4.1  were  plotted  to  determine  the 


E-plane 


H-plane 


x 
▲ 


0  =  98°  Plane 
~7 


Figure  4.1      Orientation  of  E-plane  and  H-plane 

effects  of  the  radome.  The  antenna  pattern  without  the  radome  was  plotted 
(dashed  line)  along  with  the  pattern  of  the  complete  system  (solid  line)  to 
illustrate  the  effects  of  radome  scattering. 
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1.       Code  Validation:  Scattering  from  a  Small  Disk 

The  pattern  of  an  antenna  radiating  in  the  presence  of  a  negligibly-small 
circular-disk  located  at  a  large  distance  was  computed.  Since  the  radome 
scattered  field  is  negligible,  the  predicted  output  is  that  of  a  uniformly-excited 
circular-aperture  antenna.  Figure  4.2  shows  the  plots  of  the  E-plane  pattern  for 
three  different  radius  antennas  (.5  /..  1.0  >,  and  2.0  A).  The  plots  appear  as  single 
lines  with  no  scattering.  The  gain  as  computed  by  GAIN.F,  was  compared  to  the 
theoretical  gain  of  a  circular  aperture  with  a  uniform  current  distribution 


Go- 


(2*a 


\     A      / 


(4-1) 


where  a  is  the  aperture  radius  |Ref.  10:pp.  483].   Theoretical  gains  for  the  .5  A.. 
1.0  A.  and  2.0  k  radius  apertures  are  9.94.  15.96,  and  21.98  dB  respectively.    In 
each  case  the  numerical  gain  is  within  1.5  dB  of  the  theoretical  gain. 
2.      Code  Validation:   Perfectly  Conducting  Paraboloid 

Programs  have  been  written  that  use  the  MM  solution  for  bodies  of 
revolution  developed  by  Mautz  and  Harrington  to  predict  the  patterns  and  gains 
of  reflector  antennas  [Ref.  7|.  Assuming  that  the  radome  is  a  perfectly 
conducting  parabola,  the  output  of  RADOME. F  and  GAIN.F  can  be  compared 
with  the  output  of  these  existing  programs. 

Figure  4.3  is  a  test  case  where  the  radome  is  a  parabola  with  surface 
impedance  ZL  =  8+  J0  Q.  The  dimensions  of  the  radome  matched  those  of  a 
reflector  which  was  analyzed  by  a  separate  computer  program  that  calculates 
both  pattern,  and  gain.   The  pattern  and  gain  for  the  perfectly  conducting  radome 
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cc 

73 


•20  \- 


-40- 


-60 


!    \ 

Antenna  Radius  =  .5 

Numerical  Gain  =  11.42  (dB) 
Theoretical  Gain  =  9.94  (dB) 

r 

0 


-20 


02 


-40 


-60 


0  50  LOO  150         200         250 

Theta  Angle  (deg) 


300 


Antenna  Radius  =  2.0 

Numerical  Gain  =  22.32  (dB) 
Theoretical  Gain  =  21.98  (dB) 


350 


\                          Antenna  Radius  =  1.0                              / 

\  /-^           Numerical  Gain  =  16.66  (dB)           ^-^  1 

\[      \       rheoretical  Gain  =  15.96  (dB)      /       \l 

■ 

i 

0     50    100    150    200    250    300    350 

Theta  Angle  (deg) 


0  50  100  150         200         250         300         350 

Theta  Angle  (deg) 

Figure  4.2     Scattering  from  a  Negligibly  Small  Disk 
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computed  by  RADOME.F  and  GAIN.F  are  identical  to  that  computed  for  the 
reflector. 


Feed 


a  Perfectly  Conducting  Paraboloid 


Diameter  =  5  X 


Focal  Length  =  2X 


■>-  z 


CQ 


0 
-10 
-20 
-30 
-40 
-50 
-60 


Feed  Antenna  Main 


Backscatter  from  Reflector 


Lobe 


0 


50  100 

Theta  Angle  (deg) 


.50 


Figure  4.3     Perfectly  Conducting  Paraboloid 

3.      Dielectric  Radomes 

To  determine  the  effect  of  various  radome  materials  on   antenna 
performance,  radomes  with  several  different  surface  impedances  were  analyzed. 
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Surface  impedances  were  calculated  using  (2-32)  for  a  range  of  er  and  tan<5 
typical  of  radome  materials.  The  real  and  imaginary  components  of  ZL  as  a 
function  of  er  and  tan  8  are  plotted  in  Figure  4.4.  As  shown  in  the  figure,  the 
reactance  is  essentially  independent  of  tan<S  over  the  range  8  <  tan  5  <  8.1, 
and  the  real  component  of  ZL  is  strongly  dependent  on  tan  8,  particularly  in  the 
region  of  low  er . 

A  test  radome  consisting  of  an  ogive  with  dimensions  given  in  Table  I 
and  parameters  nx  =  1/28  and  tan  8  =  8.81  was  analyzed.*    The  patterns  were 

calculated  for  dielectric  constants  ranging  2  <  er  <18,  and  the  results  are 
summarized  in  Table  III.  The  E  and  H-plane  patterns  of  an  ogive  radome  with 
Z£=34-j"1788  Q  are  shown  in  Figures  4.5  and  4.6.  The  impedance 
corresponds  to  a  radome  made  out  of  a  material  of  er  =2.8  and  tan  8  =  8.81. 
The  scattering  is  symmetric  with  respect  to  9  for  the  unscanned  antenna,  and  the 

backlobe  level  is  -45  dB  below  the  peak.  The  scanned  case  has  higher  sidelobes, 

as  expected,  with  a  lobe  of  -30  dB  at  an  angle  278°.    The  gain  loss  due  to  the 

radome  for  the  unscanned  case  is  approximately  0.1  dB,  and  for  the  scanned  case 

0.2  dB.  The  discontinuity  in  the  H-plane  patterns  at  0  =  98°  is  due  to  the 
antenna  model.    When  0  =  98°  the  field  is  (J)  polarized,  and  according  to  (2-71), 

El  has  no  COS0  factor  to  force  the  pattern  to  0  at  0  =  98°.     Note  that  the 

radome  ohmic  loss  has  not  been  included. 

Increasing  the  dielectric  constant  results  in  a  constant-shape  scattering 
pattern  with  increasing  magnitude.  Figure  4.7  is  the  E-field  pattern  for  er  =18, 
which  corresponds  to  some  ceramics.   It  is  apparent  that  the  scattering  by  the 


*  These  dimensions  are  not  intended  to  match  any  particular  existing 
radome. 
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radome  in  the  rear  hemisphere  is  significant,  and  the  gain  loss  is  nearly  4  dB. 
Table  III  shows  the  gains  with  and  without  the  radome  for  2  <  er  <  1  0.  Figure 
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Figure  4.4     Surface  Impedance  ZL  as  a  Function  of  er  and  tan  8 

4.8  is  a  plot  of  the  gain  loss  relative  to  no  radome  as  a  runction  of  er.    For 
reference  the  Fresnel  reflection  coefficient  for  a  planar  interface  is  also  shown. 
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TABLE  III.   GAIN  LOSS  FOR  AN  OGIVE  RADOME 


Z<   (Q) 

(/I    =->-              1 

A.           90  » 

Gain  (dB) 
With/Without 
Radome 

Gain  (dB) 
With/Without 
Radome 

£r 

£0 

v  tan5  =  0.01  J 

0s  =0°  U    =Q°\ 

ds  =30°  (os  =0°) 

2.0 

0-J1700 

22.21/22.32 

21.61/21.78 

2.0 

34-J1700 

22.21/22.32 

21.61/21.78 

4.0 

7.5-J560 

21.42/22.32 

20.90/21.78 

5.0 

5.3-J420 
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Figure  4.5     E-Plane  Pattern  of  Antenna  with  an  Ogive  Radome 
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Figure  4.7     E-Plane  Pattern  of  Antenna  with  an  Ogive  Radome 
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The  gain  loss  appears  to  be  linear  with  er .  This  is  probably  because  the  base  of 
the  ogive  is  open,  which  allows  reflected  waves  to  exit  the  radome.  If  a  back 
surface  was  present,  these  waves  would  be  reflected  again  causing  additional 
gain  loss. 

a.     Conical  Radome 

To  determine  the  effect  that  radome  shape  has  on  scattering,  a 
conical  radome  was  analyzed.  The  length,  width,  and  impedance  of  the  cone  was 
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Figure  4.8     Gainloss  Relative  for  an  Ogive  Radome  as  a  Function  of  er  (E- 

Plane  Scan) 


the  same  as  the  ogive  radome  discussed  previously.  A  conical  radome  presents  a 
relatively  flat  surface  that  could  possibly  cause  more  focused  scattering  than  a 
curved  shape  like  an  ogive.  Ray  tracing  of  Figure  4.9  shows  that  the  reflections 
from  the  radome  will  add  coherently  for  the  cone,  and  non-coherently  for  the 
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ogive.    However,  this  assumes  that  the  radome  is  in  the  far  field  of  the  antenna. 
The  E-plane  pattern  is  shown  in  Figure  4.10. 
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Scattered  Field 
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Figure  4.9     Scattering  from  Conical  Radome  vs.  Ogive 

A  comparison  of  Figures  4.5  and  4.10  demonstrates  that  there  is 
even  less  scattering  from  the  conical  radome  in  the  rear  hemisphere.  For  the 
scanned  case,  the  conical  radome  has  a  distinct  lobe  at  0  =  280°  due  to  the 
specular  reflection  of  the  main  lobe  from  one  side  of  the  radome.  This  is  referred 
to  as  a  reflection  lobe. 
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Figure  4.10  E-Pane  Pattern  of  Antenna  with  a  Conical 
RadomeZL  =34-  J1780  Q 
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B.    CONCLUSIONS 

This  research  has  resulted  in  two  computer  codes  that  can  be  used  to  predict 
the  pattern  degradation  and  gain  loss  due  to  an  axially  symmetric  radome  in  the 
near  field  of  a  linearly  polarized  circular  aperture.  The  computed  patterns  show 
that  the  gain  loss  is  slightly  greater  than  the  reflection  loss  for  a  planar  interface 
with  the  same  dielectric  constant  as  the  radome.  This  is  approximately  true  for 
low-loss  radomes.  For  lossy  radomes  the  impedance  has  a  very  strong  affect  on 
scattering,  while  radome  shape  has  less  effect.  Scanning  the  antenna  results  in  an 
asymmetric  scattering  pattern  and  generally  higher  sidelobe  levels. 

The  computed  gain  compares  favorably  with  values  computed  by  closed  form 
solutions  derived  from  theory,  and  the  gain  of  the  system  has  the  expected 
decrease  as  the  level  of  scattering  increases.  For  the  specific  case  of  the  ogive 
radome,  there  is  a  linear  dependence  of  the  gain  loss  to  the  relative  dielectric  of 
the  radome.  Note  that  the  gain  calculations  presented  here  do  not  include  ohmic 
loss  inside  of  the  radome. 

Improvements  in  the  original  code  include  arbitrary  aperture  illumination  and 
incorporating  mode  symmetry.  Computational  run  time  is  reduced  by  nearly  half 
when  symmetry  is  used,  and  significantly  increases  the  applicability  of  this 
method  of  analysis. 

In  order  to  apply  this  analysis  to  multilayer  radomes,  an  equivalent 
surface  impedance  must  be  determined.  The  program  can  be  modified  to  account 
for  radomes  constructed  of  nonuniform  materials  by  assigning  different 
impedances  for  each  segment  of  the  plane  curve  generating  the  BOR. 
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APPENDIX   A.    GAUSS  QUADRATURE 

Computing  the  MM  elements  of   U    and    Z    (2-19  and  2-28)  as  well  as  the 

gain  requires  the  numerical  evaluation  of  integrals.  MM  requires  numerous 
integration  in  order  to  solve  for  the  current  coefficients  /'    and  1°    of  (2-31).   Not 

*■"  "7  nj 

only  is  numerical  analysis  necessary  to  perform  the  integrals  in  the  MM,  but  also 
the  matrix  analysis  of  (2-3 1 ).  The  large  number  of  separate  integrals  necessitates 
an  integration  technique  which  provides  accuracy  in  a  few  steps  to  avoid 
excessive  run  time.  Gauss  Quadrature  is  a  method  of  numerical  integration  which 
employs  unequally  spaced  intervals  ,  and  provides  accuracy  with  a  few  points. 
[Ref.  4:pp.  154-157] 

Gauss  Quadrature  approximation  of  the  form 


l  =  ]f(H)dH  (A-l) 


by  the  sum  of  m  terms 


l  =  ^^Wkf(Hk),  (A-2) 


k  =  \ 


The  coefficients  LUk  are  weights,  and  the  Hk  s  are  abscissas  determined  from 


H=b11+b-a^  (A.3) 


where  %k'  s  are  the  m  number  of  zeros  of  the  m'th  degree  Legendre  Polynomials. 
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APPENDIX  B.  RADOME.F  LISTING 

The  material  on  the  following  pages  is  a  listing  of  the  program  RADOME.F 
described  in  chapter  III. 
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,-*****♦**********************«****«*****»*********+*************♦****»*** 


radome.f   V.3 

23  January  1992 

31  December  1992 

D.  JENN,  R.  FRANCIS,  K.  KLOPP 


C   ROGRAH 

C  DATE 

C   LATEST  REVISION 

C   PROGRAMMERS 

C 

C  NOTES: 

C    1.  ARBITRARY  CIRCULAR  APERTURE  ILLUMINATION  SPECIFIED  IN 

C      FUNCTION  TAPER. 

C    2.  NEAR  FIELD  COMPONENTS  INCLUDED.  (COMPUTED  IN  SUBROUTINE  GENEX . ) 

C    3.  ANTENNA  IS  AT  THE  BOR  COORDINATE  SYSTEM  ORIGIN. 

C   4.  LINEARLY  (X-POLARIZED)  ANTENNA. 

C    5.  CLOSED  FORM  FAR  FIELD  ANTENNA  OPTION  ADDED. 

C   6.  FLAGS  USED  THROUGHOUT:  (IMP  IS  SET  AUTOMATICALLY) 

C  IMP=0  perfect  conductor  radome  (for  test  purposes) 

C         IPRINT=0  print  pattern  points  to  unit  8 

C         ISEG=0  print  the  generating  curve  points  to  unit  8 

C         ICALC=0  compute  current  coefficients  ft  pattern 

C  #0  read  current  coefficients  from  disc  file  given 

C  by  'CURCO  and  compute  the  pattern. 

C         IFF=0  closed  form  far  field  antenna  pattern  is  used  to 

C  compute  ETF  and  EPF  (SUBROUTINE  ANTFF) 

C  #0  far  field  computed  using  CIRCRTP  at  distance  ROBS 

C         ISYM=0  use  mode  symmetry 

C  #0  turn  off  mode  symmetry  (for  testing) 

C   7.  THE  FACTOR  ETA  (=377  OHMS)  IS  OMITTED  THROUGHOUT. 

C 

C  BEGIN  MAIN  PROGRAM********************************************* 

CHARACTER  ITH 

CHARACTER*9  OUT , ETSCATTER , EPSCATTER , CURCO 

CHARACTER*12  CC 

CHARACTER*6  AN , ETFEED , EPFEED 

CHARACTER*8  GNT , GNPHI , CGP , CPOL , XPOL 

CHARACTER* 14  TPTS , PPTS , PHIPTS 

COMPLEX  EP,ET,Z( 100000) ,R(1600) ,B(800) ,C(800) ,U 

COMPLEX  UC,ET1,EP1,ET2,EP2,EC,EX,ZL0(400),ZL(2400) 

COMPLEX  CIRCR,CIRCT,CIRCP,ETF,EPF 

COMPLEX  EXP  1 , EXP2 , CON JG , CEXP , CMPLX , CT(30000 ) , IMP , JK 

DIMENSION  RH(400) ,ZH(400) ,XT(4) , AT(4) ,IPS(800) ,XP(100) ,AP(100) 

DIMENSION  A(100),X(100),EXP(500) ,ANG(500) ,ECP(500) 

DIMENSION  ECV(500) ,EXV(500) ,PHC(500) ,PHX(500) 

DIMENSION  ETSCAT(500) ,EPSCAT(500) ,ETHF(500) ,EPHF(500) 

INTEGER  CNPHI, SELECTION 

DATA  IPRINT/O/ , ICALC/1/ , IFT/O/ , ISEG/O/ , ISYM/O/ 

DATA  START, STOP/0. ,360./, PI/3. 1415926/ 

Rad=PI/180. 

ECX=0. 

BK=2.*PI 

U=(0.,1.) 

U0=(0. ,0.) 

UC=-U/4./PI 
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JK=(0. 0,6. 283185307) 
C*****  ENTER  A  LETTER  TO  INDICATE  WHICH  ITERATION  THIS  IS  *********** 
C*****    (ALL  DATA  FILES  ARE  APPENDED  WITH  THIS  LETTER)    *********** 
C*****    Cm  FILES  ARE  FOR  GENERATING  PLOTS  IN  MATLAB)    *********** 

WRITE(6,*) 'ENTER  A  LETTER  TO  INDICATE  WHICH  ITERATION  THIS  IS' 

READ(5,*)ITH 

OUT='outldbor'//ITH 

AN  ='ang'//ITH//' .m' 

ETFEED='etf ' //ITH//' .m' 

EPFEED='epf '//ITH//' .m' 

ETSCATTER= ' etscat ' //ITH// ' . m ' 

EPSCATTER='epscat'//ITH//' .m' 

CURCO='asciidat '//ITH 

CC= ' curcoef sdat ' //ITH 

CPOL='etpol'//ITH//' .m' 

XPOL='eppol'//ITH//' .m' 

IF(ICALC.EQ.O)  THEN 
WRITE(6,*) 'ENTER  THE  ANTENNA  RADIUS  (wavelengths)' 
READ (5,*)  ARAD 
C*****SELECT  THE  GEOMETRY  OF  THE  BOR************************************* 


WRITE(6,*) 
WRITE(6,*) 
WRITE(6,*) 
WRITE(6,*) 
WRITE(6,*) 
WRITE(6,*) 
WRITE(6,*) 


SELECT  BOR  GEOMETRY  BY  NUMBER.' 
NUMBER         BOR' 

1  OGIVE' 

2  SPHERE ' 

3  CONE' 

4  DISK' 

5  PARABOLA' 
READ(5,*)SELECTI0N 

IF  (SELECTION. EQ.DTHEN 

WRITE(6,*)* CALLING  SUBROUTINE  FOR  OGIVE' 
CALL  OGIVE(NP,ZH,RH,BASE,RS,ZP) 
ELSEIF(SELECTION . EQ . 2)THEN 

WRITE(6,*) 'CALLING  SUBROUTINE  FOR  SPHERE' 
CALL  TESTSPHERE ( NP , ZH , RH , BASE , RS , ZP ) 
ELSEIF(SELECTION . EQ . 3)THEN 

WRITE (6,*) 'CALLING  SUBROUTINE  FOR  CONE' 
CALL  CONE(NP,ZH(RH,BASE,RS,ZP) 
ELSEIF(SELECTION . EQ . 4)THEN 

WRITE (6,*) 'CALLING  SUBROUTINE  FOR  DISK' 
CALL  DISK(NP,ZH,RH,BASE,RS,ZP) 
ELSEIF(SELECTION . Eq . 5)THEN 

WRITE(6,*) 'CALLING  SUBROUTINE  FOR  PARABOLA' 
CALL  PARAB(NP,ZH,RH,DM,FOD(ZP) 
BASE=DM/2. 
ENDIF 
IF  (NP.GT.400)THEN 

WRITE(6,*)' MAXIMUM  NUMBER  OF  POINTS(NP)  IS  400' 
GOTO  999 
ENDIF 
ENDIF 
C************************************************************ 
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WRITE (6,*) 'ENTER  THE  FILENAMES  gaus### ' 

WRITE(6,*)  ' ' 

WRITE (6,*) 'ENTER  THE  FILENAME  IN  T  (NT)' 
READ(5,*)GNT 

WRITE (6,*) 'ENTER  THE  FILENAME  IN  PHI  (NPHI) ' 
READ(5,*)GNPHI 

WRITE(6,*) 'ENTER  THE  FILENAME  IN  X  AND  Y' 
READ(5,*)CGP 
C  OPEN  THE  FILES  FOR  THE  gaus/gaus### 
TPTS='gaus/'//GNT 
PPTS='gaus/'//GNPHI 
PHIPTS='gaus/'//CGP 
OPEN ( 1 , FILE=TPTS , STATUS= ' OLD  ' ) 
OPEN ( 2 , FILE=PPTS , STATUS= ' OLD ' ) 
0PEN(4,FILE=PHIPTS,STATUS='0LD') 
READ (1,*) NT 

IF(NT.GT.4)THEN 

WRITE (6,*) 'MAXIMUM  NUMBER  OF  POINTS(NT)  IS  4' 
GOTO  999 
ENDIF 
READ (2,*) NPHI 

IF(NPHI.GT.200)THEN 

WRITE (6,*) 'MAXIMUM  NUMBER  OF  POINTS(NPHI)  IS  200' 
GOTO  999 
ENDIF 
READ(4,*)CNPHI 

IF(CNPHI.GT.100)THEN 

WRITE(6,*) 'MAXIMUM  NUMBER  OF  POINTS(CNPHI)  IS  100' 
GOTO  999 
ENDIF 
C  LOAD  THE  WEIGHTS  AND  ABSCISSAS  IN  THE  VECTORS. 
DO  1  M=1,NT 

READ(1,*,END=1)XT(M) ,AT(M) 

1  CONTINUE 

DO  2  M=1,NPHI 

READ(2,*,END=2)X(M),A(M) 

2  CONTINUE 

DO  4  M=1,CNPHI 

READ(4,*,END=4)XP(M),AP(M) 

4    CONTINUE 

CLOSE(l) 

CL0SEC2) 

CL0SE(4) 

WRITE (6,*) 'ENTER  THE  PLOTTING  INCREMENT  IN  DEGREES' 
READ(5,*)DT 
IF(ICALC.EQ.O)  THEN 
WRITE(6,*) 'ENTER  HIGHEST  MODE' 
READ (5,*)  MODES 
NBL0CK=2*M0DES+1 
NROW=NBLOCK* (2*NP-3 ) 
C  CHECK  FOR  BOUNDS  VIOLATION 
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IFCNROW.GT. 30000)  THEN 

WRITE(6,*)  'BOUNDS  VIOLATION  FOR  ARRAY  CT( . ) ' 
GOTO  999 
ENDIF 
ENDIF 

WRITE (6,*) 'ENTER  PHI  (observation)  IN  DEGREES' 
READ(5,*)P 
PHI=P*RAD 

IF(ICALC.Eq.O)  THEN 

WRITE (6,*) 'ENTER  THE  SCAN  ANGLES  IN  DEGREES:  THETA.PHI' 
READ(5,*)  THD.PHD 
THS=THD*RAD 
PHS=PHD*RAD 

WRITE (6,*) 'ENTER  COMPLEX  IMPEDANCE,  OHMS:  (Real.Imag)' 
READ(5,*)  IMP 
ELSE 
C**********READ  CURRENT  COEFFICIENTS  IF  ICALCffO************************ 
WRITE(6,*)  'ENTER  LETTER  CODE  FOR  FILE  curcoefsdat' 
READ(5,*)  ITH 
CC= ' curcoefsdat ' //ITH 
0PEN(31,FILE=CC) 
WRITE(6,*)  'FILE  OPENED  IS  ' ,CC 
READ (31,*)  ITH , SELECTION , B ASE , NP , MT , MP , N 
READ (31,*)  NROW , MODES , NBLOCK , THS , PHS , ARAD , PHIO 
READ(31,*)  IMP.ZP.RS 
READ(31,*)  (RH(L),L=1,NP) 
READ(31,*)  (ZH(L),L=1,NP) 
READ(31,*)  (ZL0(L),L=1,NP) 
READ(31,*)  (CT(L),L=1,NR0W) 
CL0SE(31) 

WRITE(6,*)  'RUN  CODE  READ  FROM  DISC  :  ' , ITH 
ENDIF 
C**********reaDY  TO  CALCULATE  THE  PATTERN***************************** 
IT=INT((ST0P-START)/DT)+1 
MHI=M0DES+1 
0PEN(8,FILE=0UT) 

WRITE(8 , 8000 )  2 . *BASE , NP , MODES , RS , ZP 
8000  FORMAT(//, '***  BOR  RADIATION  PATTERN  FOR  A  CIRCULAR' 
*'  DISC  RADIATOR  USING  GENEX  ***', 
*//,2X,'B0R  DIAMETER  (WAVELENGTHS )=' ,F5. 2 ,/ ,2X , 

*  'NUMBER  OF  GENERATING  POINTS  (NP)= ' , 14 , / ,2X , 

*  'HIGHEST  MODE  NUMBER  USED  (MODES )= ' , 13 , 
*//,2X, 'SURFACE  RADIUS' ,F5.2,/,2X, 'ZPRIME' ,F5.2) 

WRITE(8,30)  NT.NPHI.CNPHI 
30  F0RMAT(/,12X, '  NT   NPHI   CNPHI '  , 

*  /,10X,I5,2X,I5,2X,I5) 
IF(ISEG.EQ.O)  WRITE(8,1300) 

1300   F0RMAT(/,10X, 'INDEX' ,8X, 'Z(I) ' , 10X , 'RHO(I) ' , 12X, 'SURF   IMPED') 
MP=NP-1 
MT=MP-1 
N=MT+MP 


69 


DO  52  1  =  1, NP 

IF(ABS(ZB(I)).LT. .001)  ZH(I)=0. 

IF(ABS(RH(I)).LT. .001)  RH(I)=0. 

ZHB=ZH(I)/BK 

RHB=RH(I)/BK 
C  ASSIGN  SURFACE  IMPEDANCE  AT  THIS  POINT.   THE  SURFACE  IMPEDANCE  OF  SEGMENT 
C  I  IS  ZLO(I) 

IF(ICALC.EQ.O)  Z10(I)=IMP/(120.*PI) 

IF(ISEG.EQ.O)  WRITE(8,8004)  I ,ZHB ,RHB , IMP 
52    CONTINUE 

8004   FORMAT ( 1 IX , 14 , 4X , F8 . 3 , 8X , F8 . 3 , 6X , 2F8 . 2 ) 

C************************************************************************ 
C  BIG  LOOP  *>***************************** 
C 

C  MODE  LOOP  TO  CALCULATE  THE  CURRENT  COEFFICIENTS.  POS  AND  NEG  MODES  * 
:  :ONE  IN  THE  SAME  ITERATION  OF  THE  LOOP 

C  "  v 

C  400  CONTINUE  <  **************************** 
C**********************************************ZLOAD,ZMAT,GENEX,DECOMP, SOLVE 

IF(ICALC.EQ.O)  THEN 

IF(CABSUMP).NE.O)  CALL  ZLOAD(NP  .RH.ZH  ,ZLO  ,ZL) 

DO  400  M=1,MHI 

NM=M-1 

CALL  ZMAT( NM , NM , NP , NPHI , NT , RH , ZH , X , A , XT , AT , Z ) 
C*******************************************************POSITIVE  MODE  +n 


c 
c 

t,t 

z 

t.phi   1 

z      1 

c 

+n 

'+n      .  1 

c 

ZMAT  = 

c 
c 

phi 
Z 

t 

phi, phi  1 
Z        1 

c 

+n 

+n       | 

C*********************************************************************** 
C  MODE  SYMMETRY  IN  n  FOR  ZMAT  EXISTS  BUT  IS  NOT  USED  IN  THIS  VERSION. 
C  MODE  SYMMETRY  IN  I  USED  DIRECTLY  FOR  PHIS  OF  0,90,180,  AND  270. 

C*********************************************************************** 


c 

t,t 

t  ,phi 

c 

z 

Z 

c 

-n 

-n 

c 

c 

phi  ,  t 

phi, phi 

c 

z 

Z 

c 

-n 

-n 

t.t 

t.phi 

Z 

-Z 

+n 

+n 

phi.t    phi, phi 
-Z       Z 
+n       +n 
C****************************************************CURRENT  COEFFICIENTS 

C  -1 

C 

c 
c 
c 
c 

C   CT  = 


1      ol 

0 

1      0      1 

1  I     1 

1        0 

0 

Z         1 

1    v       | 

1      -nl 

-n    1 

1      -n   1 

c 
c 
c 
c 
c 
c 

C***************************************************P0SITIVE  MODE   +n 
IF(CABSUMP)  .NE.O)    CALL  ZTOT(HT,MP  ,  ZL,  Z) 
CALL   GENEX ( NM , NP , NT , NPHI , CNPHI , XT , AT , X , A , 

*  XP.AP.THS.PHS.ARAD.RH.ZH.B) 
CALL  DECOMP(N,IPS,Z) 

CALL   SOLVE(N,IPS,Z,B,C) 

NT0P1=M0DES-NM 

NT0P2=NBL0CK-(NT0P1+1) 

NS2=NT0P1*N 

NS1=NT0P2*N 

NROW=NBLOCK*N 
C*************************************************************** 
C  STORE  CURRENT  COEFFICIENTS  IN  ONE  LONG  VECTOR  CT( . ) .  MOST  NEGATIVE 
C  MODE  IS  AT  THE  TOP;  MOST  POSITIVE  AT  THE  BOTTOM.   UNIT  30  IS  FORMATED 
C  OUTPUT;  UNIT  31  IS  FREE  FORMATED  (USED  BY  GAIN  PROGRAM). 
q* *****************************************************  ********* 

0PEN(30,FILE=CURC0) 

WRITE(30,5027)  'Mode(  '  ,NM,  ')  =  \NM,'NM  =  \NM 
DO  401  L=1,N 

WRITE(30,5026)  'CT(',L+  NS1,')  =  ' ,C(L) 

401  CT(L+NS1)=C(L) 

5027   FORMAT(/A,I,A,I/A,I/) 

5026   F0RMAT(5X,A,I4,A,F10.4,'+  (  '  ,F10 . 4, ' )*i' ) 

C***************END  OF  POSITIVE  MODE*******BEGIN  NEGATIVE  MODE  -n 
C  MODE  SYMMETRY  IS  USED  IF  ISYM=0 

IF(NM.NE.O)  THEN 

NMN=-NM 
C  USE  BRUTE  FORCE  IF  ISYM  #  0 

IF(ISYM.NE.O)  THEN 

CALL  GENEX (NMN , NP , NT , NPHI , CNPHI , XT , AT , X , A , 

*  XP,AP,THS(PHS(ARAD,RH,ZH,B) 

CALL  ZMAT( NMN , NMN , NP , NPHI , NT , RH , ZH , X , A , XT , AT , Z ) 

IF(CABSUMP).NE.O)  CALL  ZTOT(MT,MP  ,ZL,Z) 

CALL  DECOMP(N,IPS,Z) 

CALL  SOLVE(N,IPS,Z,B,C) 

ENDIF 

WRITE(30,5027),Mode(\NMN,')  = 

DO  402  L=1,MT 

CT(L+NS2)=+C(L) 

WRITE(30,5026)  'CT( ' .L+NS2, ' )  = 

402  CONTINUE 
DO  403  L=MT+1,N 
IF(ISYM.NE.O)  CT(L+NS2)=+C(L) 
IF(ISYM.Eq.O)  CT(L+NS2)=-C(L) 
WRITE(30,5026)  'CT( ' .L+NS2 , ' )  =  \CT(L+NS2) 


, NMN , ' NM  =  ' , NM 


,CT(L+NS2) 


403     CONTINUE 

C******************************************END  OF  NEGATIVE  MODE  -n 

ENDIF 
400   CONTINUE 

CL0SEC30) 
(2**+***********************************+******************+*********+*** 

C  WRITE  THE  VECTOR  OF  CURRENT  COEFFICIENTS  ON  DISC  FOR  GAIN  CALC. 

C  (FREE  FORMAT  TO  UNIT  31.) 

C*********************************************************************** 

0PEN(31,FILE=CC) 

WRITE (31,*)  ITH , SELECTION , BASE , NP , MT , MP , N 

WRITE (31,*)  NROW , MODES , NBLOCK , THS , PHS , ARAD , PHI 

WRITE(31,*)  IMP.ZP.RS.IFF 

WRITE(31,*)  (RH(L),L=1,NP) 

WRITE(31,*)  (ZH(L),L=1,NP) 

WRITE(31,*)  (ZLO(L) ,L=1,NP) 

WRITE(31,*)  (CT(L),L=1,NR0W) 

WRITE(31,*)    GNT.GNPHI.CGP 

WRITE(6,*)    'CURRENT  COEFFICIENTS   WRITTEN  ON  DISC 

CL0SE(31) 

ENDIF 
C*********************************** ************************ ************* 
C  END  OF  BIG  LOOP   *>************************* 

C  *  v 

C*******<   ************   *************** 

c************************************************************************ 
c 

C  BEGIN  PATTERN  LOOP  FROM  START  TO  STOP  IN  INCREMENTS  OF  DT  (ALL  IN  DEG) 
C 

DO  500  1  =  1, IT 
THETA=START  +  DT*FL0AT(I-1 ) 
THX=THETA*RAD 
PHR=PHI 

IF(THETA.GT.180.)  THEN 

PHR=PHI+PI 

THX=(360.-THETA)*RAD 

ENDIF 

ET1=(0. ,0.) 

EP1=(0. ,0.) 

ET2=(0. ,0. ) 

EP2=(0. ,0.) 

DO  300  M=1,MHI 

NM=M-1 

EXP1=CEXP(CMPLX(0. ,FLOAT(NM)*PHR) ) 

EXP2=C0NJG(EXP1) 
C***********************************************************PLANE 

CALL  PLANE ( NM , NM , NP , NT , RH , ZH , XT , AT , THX , R ) 

NT0P1=M0DES-NM 

NT0P2=NBL0CK-(NT0P1+1) 

NS2=NT0P1*N 

NS1=NT0P2*N 
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DO  250  L=1,MT 

ET1=ET1+R(L)*CT(L+NS1)*EXP1 

EP1=EP1+R(L+N)*CT(L+NS1)*EXP1 

IF(NM.EQ.O)  GOTO  250 

ET1=ET1+R(L)*CT(L+NS2)*EXP2 

EP1=EP1-R(L+N)*CT(L+NS2)*EXP2 
250   CONTINUE 

DO  260  L=1,MP 

ET2=ET2+R(L+MT)*CT(L+NS1+MT)*EXP1 

EP2=EP2+R(L+MT+N)*CT(L+NS1+MT)*EXP1 

IF(NM.EQ.O)    GO   TO    260 

ET2=ET2-R(L+MT)*CT(L+NS2+MT)*EXP2 

EP2=EP2+R(L+MT+N)*CT(L+NS2+MT)*EXP2 
260   CONTINUE 
300   CONTINUE 

C  FEED  CONTRIBUTION  IN  THE  FAR  FIELD  IS  ETF.EPF 
C******** ******************************************* *BESS J 1 . 

C  USE  CLOSED  FORM  FEED  EXPRESSION  FROM  SUBROUTINE  ANTFF  IF 
C  IFF=0;  OTHERWISE  USE  BRUTE  FORCE  EVALUATION  FROM  CIRCRTP 
C 

R0BS=1000. 
IF(IFF.EQ.O)  THEN 

CALL  ANTFF (THX , PHR , THS , PHS , ARAD , ETF , EPF ) 
ELSE 
RHB=ROBS*SIN(THX) 
ZHB=ROBS*COS(THX) 

CALL  CIRCRTP (CNPHI , XP , AP , ARAD , THS , PHS , 
*  _    PHR.RHB.ZHB.CIRCR.CIRCT.CIRCP) 

C  REMOVE  THE  1/R  DEPENDENCE  BECAUSE  EXP(-jkR)/R  IS  OMITTED  IN 
C  THE  SCATTERED  FIELDS  ET  AND  EP 
ETF=CIRCT*ROBS 
EPF=CIRCP*ROBS 
ENDIF 
C**************************************************** 
ETHF(I)=CABS(ETF) 
EPHF(I)=CABS(EPF) 
ET=(ET1+ET2)*UC 
EP=(EP1+EP2)*UC 
ETSCAT(I)=CABS(ET) 
EPSCAT(I)=CABS(EP) 
C  TOTAL  E-THETA  AND  E-PHI  COMPONENTS 
EC=ET+ETF 
EX=EP+EPF 
ECV(I)=CABS(EC) 
EXV(I)=CABS(EX) 
ECR=REAL(EC) 
ECI=AIMAG(EC) 
EXR=REAL(EX) 
EXI=AIMAG(EX) 

PHC(I)=ATAN2(ECI,ECR+l.e-10)/RAD 
PHX(I)=ATAN2(EXI,EXR+l.e-10)/RAD 


ANG(I)=THETA 

ECX=AMAX1(ECX,ECV(I) ,EXV(I)) 
500  CONTINUE 

WRITE (6,*)  'MAX  E  VALUE=',ECX 

WRITE (8 , 103 )  P , THS/RAD , PHS/RAD , ARAD , ECX 
103  F0RMAT(/,10X, 'PHI  OF  RECEIVER  (DEG)= ' ,F8 . 2,/ , 10X , 

*  'ANTENNA  SCAN:  THETA  (DEG)  =  '  ,F8 . 2 ,/ , 10X, 

*  '  PHI    (DEG)  =  \F8.2,/,10X, 

*  'ANTENNA  RADIUS:       ARAD= ' ,F8 . 3 , / , 10X , 

*  'MAXIMUM  FIELD  VALUE  (V/M)=  '  ,E15 . 5) 
IF(IFF.EQ.O)  THEN 

WRITE(8,113) 

ELSE 
WRITE(8,213)  ROBS 

ENDIF 
113  F0RMAT(/,10X, 'CLOSED  FORM  FAR  FIELD  PATTERN  FROM  ANTFF  IS  USED') 
213  F0RMAT(/,10X, 'FAR  FIELD  COMPUTED  USING  CIRCRTP,  ROBS= ' ,F10. 2) 
C  STORE  DATA  FOR  NORMALIZED  CO-  AND  CROS-POLARIZED  PATTERNS 

DO  600  1=1, IT 

ECP(I)=(ECV(I)/ECX)**2 

EXP(I)=(EXV(I)/ECX)**2 

ECP(I)=AMAX1(ECP(I) , l.E-6) 

EXP(I)=AMAX1(EXP(I) , l.E-6) 

ECP(I)=10.*AL0G10(ECP(I)) 

EXP(I)  =  10.*AL0G10(EXP(D) 
600  CONTINUE 

IF(IPRINT.Eq.O)  THEN 

WRITE(8,5015) 

5015  FORMAT (//,7X, 'ANGLE' ,15X, 'CO-POL ' , 25X , 'X-POL' ,/,7X, 
l'(DEG) ' ,4X,2(* (VOLTS)' ,4X, '(DEG) ' ,3X, ' (DB-REL) ' ,4X)) 

DO  9000  L=1,IT 

WRITE(8,5016)  ANG(L) ,ECV(L) ,PHC(L) ,ECP(L) ,EXV(L) ,PHX(L) 
1,EXP(L) 

5016  F0RMAT(5X,F6.2(3X,2(F8.4,3X,F7.2,3X)F7.2,3X)) 
9000  CONTINUE 

ENDIF 

0PEN(2,file=AN) 
0PEN(3,file=CP0L) 
0PEN(4,file=XP0L) 
0PEN(7,file=ETSCATTER) 
0PEN(8,file=EPSCATTER) 
0PEN(9,file=ETFEED) 
0PEH(10,file=EPFEED) 
DO  9097  1=1, IT 
WRITE(2, 5019)1,  ANG(I) 
WRITEC3, 5020)1,  ECP(I) 
WRITE(4, 5021)1,  EXP(I) 
WRITE(7, 5022)1,  ETSCAT(I) 
WRITE (8,5023) I,  EPSCAT(I) 
WRITE(9, 5024)1,  ETHF(I) 
9097  WRITE(10,5025)I,  EPHF(I) 
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5019  FORMAT (  '  ANG(  '  ,  I ,  '  )=  \F8.  3  ,  '  ;  ' ) 

5020  FORMAT ( ' ECP(  ',1 ,  ' )= ' ,F8. 3 ,  ' ; ' ) 

5021  FORMATCEXPC  ,  I ,  '  )=  '  ,F8  .  3  ,  '  ;  ' ) 

5022  FORMAT(  'ETSCATC  ,1,  '  )  =  '  ,F8.3,  '  ;  '  ) 

5023  FORMATCEPSCATC  ,  I ,  '  )  =  '  ,F8  .  3,  '  ;  •  ) 

5024  FORMAT (  'ETHFC  ,  I ,  '  )  =  '  ,F8  .  3  ,  '  ;  '  ) 

5025  FORMATCEPHFC  ,  I ,  '  )=  '  ,F8  .  3  ,  '  ;  ' ) 
CL0SE(2) 

CL0SE(3) 

CL0SEC4) 

CL0SE(7) 

CL0SE(8) 

CLOSEO) 

CLOSE(IO) 
999  CONTINUE 

STOP 

END 
C********************************************************END  OF  MAIN  PROGRAM. 
C 

C********************************************************SUBROUTINE  ZMAT. 
C  REFERENCE:  AN  IMPROVED  E-FIELD  SOLUTION  FOR  CONDUCTING  BOR 
C  J.R.MAUTZ  AND  R.F.  HARRINGTON 

C  TECHNICAL  REPORT  TR-80-1 

C  ROME  AIR  DEVELOPMENT  CENTER 

C  CONTRACT  NO.  F-30602-79-C-0011 

SUBROUTINE  ZMAT ( M 1 , M2 , NP , NPHI , NT , RH , ZH , X , A , XT , AT , Z ) 


C  COMPUTE  THE  MM  IMPEDANCE  MATRIX  ELEMENTS.   THIS  IS  FROM  MAUTZ  AND 

C  HARRINGTON  (NO  CHANGES). 

C 

COMPLEX  Z( 100000) ,U1 ,U2,U3 ,U4 ,U5 ,U6 ,U7 ,U8 ,U9 ,UA,UB ,G4A(4) , G5A(4) 

COMPLEX  G6A(4) ,G4B(4) ,G5B(4) ,G6B(4) ,H4A,H5A,H6A,H4B ,H5B 

COMPLEX  H6B(UC,UD,GA(100),GB(100) 

DIMENSION  RH(400) ,ZH(400) ,X(100) , A( 100) ,XT(4) , AT(4) 

DIMENSION  RS(400) ,ZS(400) ,D(400) ,DR(400) ,DZ(400) 

DIMENSION  DM(400) ,C2( 100) ,C3( 100) ,R2(4) ,Z2(4) 

DIMENSION  C4(100) ,C5(100) ,C6(100) ,Z7(4) ,R7(4) ,Z8(4) ,R8(4) 

CT=2. 

CP=.l 

DO  10  1=2, HP 

12=1-1 

RS(I2)=.5*(RH(I)+RH(I2)) 

ZS(I2)=.5*(ZH(I)+ZH(I2)) 

D1=.5*(RH(I)-RH(I2)) 

D2=.5*(ZH(I)-ZH(I2)) 

D(I2)=SqRT(Dl*Dl+D2*D2) 

DR(I2)=D1 

DZ(I2)=D2 

IF(RS(I2).Eq.0.)  RS(I2)=1. 

DM(I2)=D(I2)/RS(I2) 
10  CONTINUE 


M3=M2-M1+1 
M4=M1-1 
PI2=1. 570796 

DO  11  K=1,NPHI 

PH=PI2*(X(K)+1.  ) 

C2(K)=PH*PH 

SN=SIN(.5*PH) 

C3(K)=4.*SN*SN 

A1=PI2*A(K) 

D4=.5*A1*C3(K) 

D5=A1*C0S(PH) 

D6=A1*SIN(PH) 

M5=K 

DO  29  M=1,M3 

PHM=(M4+M)*PH 

A2=C0S(PHM) 

C4(H5)=D4*A2 

C5(H5)=D5*A2 

C6(M5)=D6*SIN(PHM) 

M5=M5+NPHI 
29  CONTINUE 
11  CONTINUE 

MP=NP-1 

MT=MP-1 

N=MT+MP 

N2N=MT*N 

N2=N*N 

Ul=(0.,.5) 

U2=(0. ,2.) 

JN=-1-N 

DO  15  JQ=1,MP 

Kq=2 

IF(JQ.EQ.l)  KQ=1 

IF(JQ.EQ.MP)  KQ=3 

R1=RS(JQ) 

Z1=ZS(JQ) 

D1=D(JQ) 

D2=DR(jq) 

D3=DZ(jq) 

D4=D2/R1 

D5=DM(jq) 

SV=D2/D1 

CV=D3/D1 

T6=CT*D1 

T62=T6+D1 

T62=T62*T62 

R6=CP*R1 

R62=R6*R6 

DO  12  L=1,NT 

R2(L)=R1+D2*XT(L) 

Z2(L)=Z1+D3*XT(L) 
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12  CONTINUE 
U3=D2*U1 
U4=D3*U1 
DO  16  IP=1,MP 
R3=RS(IP) 
Z3=ZS(IP) 
R4=R1-R3 
Z4=Z1-Z3 
FM=R4*SV+Z4*CV 
PHM=ABS(FM) 
PH=ABS(R4*CV-Z4*SV) 
D6=PH 

IF(PHM.LE.Dl)  GO  TO  26 
D6=PHM-D1 
D6=SQRT(D6*D6+PH*PH) 

26  IF(IP.EQ.JQ)  GO  TO  27 
KP=1 

IF(T6.GT.D6)  KP=2 
IF(R6.GT.D6)  KP=3 

GO  TO  28 

27  KP=4 

28  GO  TO  (41,42,41,42) ,KP 
42  DO  40  L=1,NT 

D7=R2(L)-R3 

D8=Z2(L)-Z3 

Z7(L)=D7*D7+D8*D8 

R7(L)=R3*R2(L) 

Z8(L)=.25*Z7(L) 

R8(L)=.25*R7(L) 
40  CONTINUE 

Z4=R4*R4+Z4*Z4 

R4=R3*R1 

R5=.5*R3*SV 

DO  33  K=1,NPHI 

A1=C3(K) 

RR=Z4+R4*A1 

UA=0. 

UB=0. 

IF(RR.LT.T62)  GO  TO  34 

DO  35  L=1,NT 

R=SqRT(Z7(L)+R7(L)*Al) 

SN=-SIN(R) 

CS=COS(R) 

UC=AT(L)/R*CMPLX(CS,SN) 

UA=UA+UC 

UB=XT(L)*UC+UB 
35  CONTINUE 

GO  TO  36 
34  DO  37  L=1,NT 

R=SQRT(Z8(L)+R8(L)*A1) 

SN=-SIN(R) 
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CS=COS(R) 

UC=AT(L)/R*SN*CMPLX(-SN,CS) 

UA=UA+UC 

UB=XT(L)*UC+UB 

37  CONTINUE 
A2=FM+R5*A1 
D9=RR-A2*A2 
R=ABS(A2) 
D7=R-D1 
D8=R+D1 

D6=SQRT(D8*D8+D9) 
R=SQRT(D7*D7+D9) 
IF(D7.GE.O. )  GO  TO  38 
A1=AL0G((D8+D6)*(-D7+R)/D9)/D1 
GO  TO  39 

38  A1=AL0G((D8+D6)/(D7+R))/D1 

39  UA=A1+UA 
UB=A2*(4./(D6+R)-A1)/D1+UB 

36  GA(K)=UA 

GB(K)=UB 
33  CONTINUE 

K1=0 

DO  45  M=1,M3 

H4A=0. 

H5A=0. 

H6A=0. 

B4B=0. 

H5B=0. 

H6B=0. 

DO  46  K=1,NPHI 

K1=K1+1 

D6=C4(K1) 

D7=C5(K1) 

D8=C6(K1) 

UA=GA(K) 

UB=GB(K) 

H4A=D6*UA+B4A 

H5A=D7*UA+H5A 

H6A=D8*UA+H6A 

B4B=D6*UB+H4B 

H5B=D7*UB+B5B 

H6B=D8*UB+H6B 
46  CONTINUE 

G4A(M)=B4A 

G5A(M)=B5A 

G6A(M)=B6A 

G4B(M)=B4B 

G5B(M)=B5B 

G6B(M)=H6B 
45  CONTINUE 

IF(KP.NE.4)  GO  TO  47 
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A2=D1/(PI2*R1) 

D6=2./D1 

D8=0. 

DO    63   K=1,NPHI 

A1=R4*C2(K) 

R=R4*C3(K) 

IF(R.LT.T62)    GO  TO   64 

D7=0. 

DO   65   L=1,NT 

D7=D7+AT(L)/SqRT(Z7(L)+Al) 

65  CONTINUE 
GO  TO  66 

64  A1=A2/(X(K)+1.) 

D7=D6*AL0G(A1+SQRT(1 . +A1*A1) ) 

66  D8=D8+A(K)*D7 
63  CONTINUE 

A1=.5*A2 

A2=1./A1 

D8=-PI2*D8+2./Rl*(BL0G(A2)+A2*BL0G(Al)) 

DO  67  M=1,H3 

G5A(M)=D8+G5A(M) 

67  CONTINUE 
GO  TO  47 

41  DO  25  M=1,M3 

G4A(M)=0. 

G5A(M)=0. 

G6A(M)=0. 

G4B(M)=0. 

G5B(M)=0. 

G6B(M)=0. 
25  CONTINUE 

DO  13  L=1,NT 

A1=R2(L) 

R4=A1-R3 

Z4=Z2(L)-Z3 

Z4=R4*R4+Z4*Z4 

R4=R3*A1 

DO  17  K=1,NPHI 

R=SQRT(Z4+R4*C3(K)) 

SN=-SIN(R) 

CS=COS(R) 

GA(K)=CMPLX(CS,SN)/R 
17  CONTINUE 

D6=0. 

IF(R62.LE.Z4)  GO  TO  51 

DO  62  K=1,NPHI 

D6=D6+A(K)/SQRT(Z4+R4*C2(K)) 
62  CONTINUE 

Z4=3 . 141593/SQRT(Z4/R4) 

D6=-PI2*D6+AL0G(Z4+SQRT( 1 . +Z4*Z4) )/SQRT(R4) 
51  A1=AT(L) 
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A2=XT(L)*A1 

K1=0 

DO  30  M=1,M3 

U5=0. 

U6=0. 

U7=0. 

DO  32  K=1,NPHI 

UA=GA(K) 

K1=K1+1 

U5=C4(K1)*UA+U5 

U6=C5(K1)*UA+U6 

U7=C6(K1)*UA+U7 
32  CONTINUE 

U6=D6+U6 

G4A(M)=A1*U5+G4A(M) 

G5A(M)=A1*U6+G5A(M) 

G6A(M)=A1*U7+G6A(M) 

G4B(M)=A2*U5+G4B(M) 

G5B(M)=A2*U6+G5B(M) 

G6B(M)=A2*U7+G6B(M) 
30  CONTINUE 
13  CONTINUE 
47  A1=DR(IP) 

UA=A1*U3 

UB=DZ(IP)*U4 

A2=D(IP) 

D6=-A2*D2 

D7=D1*A1 

D8=D1*A2 

JM=JN 

DO  31  M=1,M3 

FM=M4+M 

A1=FM*DM(IP) 

H5A=G5A(M) 

H5B=G5B(H) 

H4A=G4A(M)+B5A 

H4B=G4B(M)+H5B 

H6A=G6A(M) 

H6B=G6B(M) 

U7=UA*H5A+UB*H4A 

U8=UA*H5B+UB*H4B 

U5=U7-U8 

U6=U7+U8 

U7=-U1*H4A 

U8=D6*H6A 

U9=D6*H6B-A1*H4A 

UC=D7*(B6A+D4*H6B) 

UD=FH*D5*H4A 

K1=IP+JM 

K2=K1+1 

K3=K1+N 
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K4=K2+N 

K5=K2+MT 

K6=K4+MT 

K7=K3+N2N 

K8=K4+N2N 

GO  TO  (18, 20, 19), KQ 

18  Z(K6)=U8+U9 
IF(IP.EQ. 1)  GO  TO  21 
Z(K3)=Z(K3)+U6-U7 
Z(K7)=Z(K7)+UC-UD 
IF(IP.EQ.MP)  GO  TO  22 

21  Z(K4)=U6+U7 
Z(K8)=UC+UD 
GO  TO  22 

19  Z(K5)=Z(K5)+U8-U9 
IF(IP.EQ. 1)  GO  TO  23 
Z(K1)=Z(K1)+U5+U7 
Z(K7)=Z(K7)+UC-UD 
IF(IP.EQ.HP)  GO  TO  22 

23  Z(K2)=Z(K2)+U5-U7 
Z(K8)=UC+UD 

GO  TO  22 

20  Z(K5)=Z(K5)+U8-U9 
Z(K6)=U8+U9 
IF(IP.EQ.l)  GO  TO  24 
Z(K1)=Z(K1)+U5+U7 
Z(K3)=Z(K3)+U6-U7 
Z(K7)=Z(K7)+UC-UD 
IF(IP.EQ.MP)  GO  TO  22 

24  Z(K2)=Z(K2)+U5-U7 
Z(K4)=U6+U7 
Z(K8)=UC+UD 

22  Z(K8+MT)=U2*(D8*(H5A+D4*H5B)-A1*UD) 
JM=JM+N2 

31  CONTINUE 
16  CONTINUE 

JN=JN+N 
15  CONTINUE 

RETURN 

END 
C* ******************************************* ************SUBROUTINE  SOLVE. 
C  REFERENCE  :  TECHNICAL  REPORT  TR-80-1 

SUBROUTINE  SOLVE(N, IPS ,UL,B,X) 
C 

C  SEE  MAUTZ  k   HARRINGTON  FOR  DETAILS 
C 

COMPLEX  UL(100000),B(800),X(800),SUM 

DIMENSION  IPS(800) 

NP=N+1 

IP=IPS(1) 

X(1)=B(IP) 
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DO  2  1=2, H 
IP=IPS(I) 
IPB=IP 
IM1=I-1 

SUM=(0.  ,0.) 

DO  1  J=1,IM1 

SUM=SUM+UL(IP)*X(J) 

1  IP=IP+N 

2  X(I)=B(IPB)-SUM 
K2=N*(N-1) 
IP=IPS(N)+K2 
X(N)=X(N)/UL(IP) 
DO  4  IBACK=2,N 
I=NP-IBACK 
K2=K2-N 
IPI=IPS(I)+K2 
IP1=I+1 
SUH=(0. ,0.) 
IP=IPI 

DO  3  J=IP1,N 
IP=IP+N 

3  SUM=SUM+UL(IP)*X(J) 

4  X(I)=(X(I)-SUM)/UL(IPI) 
RETURN 

END 
C ************************* *********************SUBROUTINE  DECOMP 
C  REFERENCE:  TECHNICAL  REPORT  TR-80-1 

SUBROUTINE  DECOMP(N , IPS.UL) 
C 

C  SEE  MAUTZ  4  HARRINGTON  FOR  DETAILS 
C 

COMPLEX  UL( 100000) .PIVOT, EM 

DIMENSION  SCL(800),IPS(800) 

DO  5  1*1,1 

IPS(I)=I 

RN=0. 

J1=I 

DO  2  J=1,N 

ULM=ABS(REAL(UL(J1)))+ABS(AIMAG(UL(J1))) 

J1=J1+N 

IF(RN-ULM)  1,2,2 

1  RN=ULM 

2  CONTINUE 
SCL(I)=1./RN 

5  CONTINUE 
NM1=N-1 
K2=0 

DO  17  K=1,NM1 

BIG=0. 

DO  11  I=K,N 

IP=IPS(I) 
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IPK=IP+K2 

SIZE=(ABS(REAL(UL(IPK)))+ABS(AIMAG(UL(IPK))))*SCL(IP) 
IF(SIZE-BIG)  11,11,10 

10  BIG=SIZE 
IPV=I 

11  CONTINUE 
IF(IPV-K)  14,15,14 

14  J=IPS(K) 
IPS(K)=IPS(IPV) 
IPS(IPV)=J 

15  KPP=IPS(K)+K2 
PIVOT=UL(KPP) 
KP1=K+1 

DO  16  I=KP1,N 
KP=KPP 

IP=IPS(I)+K2 
EM=-UL(IP)/PIVOT 
18  UL(IP)=-EM 
DO  16  J=KP1,N 
IP=IP+N 
KP=KP+N 
UL(IP)=UL(IP)+EM*UL(KP) 

16  CONTINUE 
K2=K2+N 

17  CONTINUE 
RETURN 
END 

C*************************************************FUNCTION  BLOG 
C  REFERENCE:  TECHNICAL  REPORT  TR-80-1  ; (page  56) 

FUNCTION  BLOG(X) 

IFCX.GT. . 1)  GO  TO  1 

X2=X*X 

BLOG=(( .075*X2-. 1666667)*X2+1. )*X 

RETURN 
1  BL0G=AL0G(X+SQRT(1.+X*X)) 

RETURN 

END 
C**** **************************************** ********SUBROUTINE  PLANE 
C  REFERENCE:  TECHNICAL  REPORT  TR-80-1  ; (pages  57-64) 

SUBROUTINE  PLANE (Ml , M2 , NP , NT , RH , ZH , XT , AT , THR , R) 
C 

C  PLANE  WAVE  EXCITATION  VECTOR  IN  THE  FAR  FIELD.  FROM  MAUTZ  AND  HARRINGTON. 
C  MODIFIED  TO  DO  ONLY  ONE  ANGLE  AND  FREQUENCY  PER  CALL. 
C  **  NEW  BESSEL  FUNCTION  ROUTINE  FROM  NUM.  RECIPES  ** 
C 

COMPLEX  R(1600),U,U1,UA,UB,FA(1500),FB(1500),F2A,F2B,F1A,F1B 

COMPLEX  U2.U3.U4.U5 

DIMENSION  RH(400) ,ZH(400) ,XT(4) , AT(4) ,R2(4) ,Z2(4) 

MP=NP-1 

MT=MP-1 

N=MT+MP 
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N2=2*N 

CC=COS(THR) 

SS=SIN(THR) 

U=(0. ,1.) 

Ul=3. 141593*U**M1 

M3=M1+1 

M4=M2+3 

IF(Ml.EQ.O)  M3=2 

M5=Ml+2 

M6=M2+2 

DO  12  IP=1,MP 

K2  =  IP 

I=IP+1 

DR=.5*(RH(I)-RH(IP)) 

DZ=.5*(ZH(I)-ZH(IP)) 

D1=SQRT(DR*DR+DZ*DZ) 

R1=.25*(RH(I)+RH(IP)) 

IF(Rl.EQ.O.)  Rl=l. 

Z1=.5*(ZH(I)+ZH(IP)) 

DR=.5*DR 

D2=DR/R1 

DO  13  L=1,MT 

R2(L)=R1+DR*XT(L) 

Z2(L)=Z1+DZ*XT(L) 
13  CONTINUE 

D3=DR*CC 

D4=-DZ*SS 

D5=D1*CC 

DO  23  M=M3,M4 

FA(M)=0. 

FB(M)=0. 
23  CONTINUE 

DO  15  L=1,NT 

X=SS*R2(L)*2. 

ARG=Z2(L)*CC 

UA=AT(L)*CMPLX(COS(ARG),SIN(ARG)) 

UB=XT(L)*UA 
C  THIS  LINE  REPLACES  HARRINGTON'S  BLOCK 

DO  25  M=M3,H4 

BES=BESSJ(H-2,X) 

FA(M)=BES*UA+FA(M) 

FB(M)=BES*UB+FB(H) 

25  CONTINUE 
15  CONTINUE 

IF(Ml.NE.O)  GO  TO  26 

FA(1)=-FA(3) 

FB(1)=-FB(3) 

26  UA=U1 

DO  27  M=M5,M6 

M7=M-1 

M8=M+1 


84 


F2A=UA*(FA(M8)+FA(M7)) 

F2B=UA*(FB(M8)+FB(M7)) 

UB=U*UA 

F1A=UB*(FA(M8)-FA(M7)) 

F1B=UB*(FB(M8)-FB(M7)) 

U4=D4*UA 

U2=D3*F1A+U4*FA(M) 

U3=D3*F1B+U4*FB(M) 

U4=DR*F2A 

U5=DR*F2B 

K1=K2-1 

K4=K1+N 

K5=K2+N 

R(K2+MT)=-D5*(F2A+D2*F2B) 

R(K5+MT)=D1*(F1A+D2*F1B) 

IF(IP.EQ.l)  GO  TO  21 

R(K1)=R(K1)+U2-U3 

R(K4)=R(K4)+U4-U5 

IF(IP.Eq.HP)  GO  TO  22 

21  R(K2)=U2+U3 
R(K5)=U4+U5 

22  K2=K2+N2 
UA=UB 

27  CONTINUE 
12  CONTINUE 

RETURN 

END 
C******** *************************************** *****SUBR0UTINE  ZLOAD , 
C  REFERENCE:  COMPUTATION  OF  RADIATION  AND  SCATTERING  FOR 
C  LOADED  BODIES  OF  REVOLUTION 

C  HARRINGTON  AND  MAUTZ 

C  REPORT  AFCRL-70-0046 

C  AIRFORCE  CAMBRIDGE  RESEARCH  LABS 

C  CONTRACT  NO.  F-19628-68-C-0180 

SUBROUTINE  ZLOAD(NP ,RH ,ZH ,ZO,Z) 
C 

C  COMPUTES  IMPEDANCE  MATRIX  ELEMENTS  FOR  LOADED  BODIES  OF  REV 
C  ZO(I)  IS  THE  SURF  IMPEDANCE  OF  THE  ITH  SEGMENT  (NP-1  SEGMENTS) 
C  Z(.)  ARE  THE  IMPEDANCE  MATRIX  TERMS  (TRIDIAGONAL  FOR  T-T 
C  SUBMATRIX;  DIAGONAL  FOR  P-P  SUBMATRIX).  STORED  IN  COL  VECTOR. 
C 

COMPLEX  C1,C2,ZO(400),Z(2400),X1,X2,X3,Y1,Y2,Y3,FN(400) 

COMPLEX  U1.U2.U3, XI, YI 

DIMENSION  RH(400) ,ZH(400) ,RS(400) ,D(400) ,SV(400) 

PI=3. 14159 

MT=NP-2 

MP=NP-1 

N=MT+MP 

DO  10  IP=2,NP 

II=IP-1 

DR=RH(IP)-RH(II) 
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DZ=ZH(IP)-ZH(II) 

D(II)=SQRT(DR*DR+DZ*DZ) 
SV(II)=DR/D(II) 
RS(II)  =  .5*(RH(IP)+RH(ID) 
DS=D(II)*SV(II)/2. 
Q1=RS(II)+DS 
Q2=RS(II)-DS 
FN(II)=1. 

IF((ABS(Q2) .GT.1.E-6).AND. (ABS(Ql) .GT.l.E-6)) 
*  FH(II)=AL0G(Q1/Q2) 
10  CONTINUE 
L0=MT*3-2 
DO  20  1=1, MP 
C1=PI*Z0(I) 
IF(I.EQ.MP)  GO  TO  80 
KI=2 

IF(I.EQ.l)  KI=1 
IF(I.EQ.MT)  KI=3 
11=1+1 

C2=PI*Z0(II) 
A=SV(I) 

IF(ABS(A) .LT.l.E-6)  GO  TO  41 
Xl=Cl*FH(I)/2./A 

X2=C1*2./A*(1.-RS(I)*FN(I)/D(I)/A) 
X3=-X2*RS(I)/D(I)/A 
GO  TO  42 

41  CONTINUE 
Xl=Cl/2./RS(I)*D(I) 
X2=(0. ,0.) 
X3=Cl*D(I)/6./RS(I) 

42  CONTINUE 
A=SV(II) 

IF(ABS(A). LT.l.E-6)  GO  TO  45 

Yl=C2*FN(II)/2./A 

Y2=C2*2./A*(1.-RS(II)*FN(II)/D(II)/A) 

Y3=-Y2*RS(II)/D(II)/A 

GO  TO  40 
45  CONTINUE 

Yl=C2/2./RS(II)*D(II) 

Y2=(0. ,0.) 

Y3=C2*D(II)/6./RS(II) 
40  CONTINUE 
C 

C  DEFINE  TRIDIAGONAL  ELEMENTS  FOR  T-T  SUBMATRIX  (STORED  IN  COLS) 
C  (Ul-  DIAG;  U2-  LOWER;  U3-  UPPER) 
C 

XI=X1+X2+X3 

YI=Y1-Y2+Y3 

IF(RH(I). LT.l.E-6)  XI=C1/SV(I) 

IF(RH(II). LT.l.E-6)  YI=C2/SV(II) 

U1=XI+YI 
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U2=X1-X3 

U3=Y1-Y3 

L=2+(I-2)*3 

IF(KI.Eq.l)  L=0 

L1=L+1 

L2=L+2 

L3=L+3 

go  to  (50,60,70) ,ki 
50  Z(L1)=U1 

Z(L2)=U2 

GO  TO  80 
60  Z(L1)=U3 

Z(L2)=U1 

Z(L3)=U2 

GO  TO  80 
70  Z(L1)=U3 

Z(L2)=U1 
80  Z(L0+I)=2.*C1*D(I)/RS(I) 
20  CONTINUE 

RETURN 

END 

SUBROUTINE  ZT0T(MT,MP ,ZL,Z) 
C 

C  ADDS  THE  SURF  IMPEDANCE  TERMS  TO  THE  TRIDIAGONAL  ELEMENTS  OF 
C  THE  BOR  IMPEDANCE  MATRIX  Z. 
C 

COMPLEX  ZL(2400) .Z(IOOOOO) 

N=MT+MP 

M0=MT*3-2 

DO  100  1=1, MP 

L0=MT*N+(I-1)*N+MT 

IF(I.EQ.MP)  GO  TO  80 

KI=2 

IF(I.Eq.l)  KI=1 

IF(I.EQ.MT)  KI=3 

L2=(I-1)*N+I 

L1=L2-1 

L3=L2+1 

M=2+3*(I-2) 

IF(KI.EQ.l)  M=0 

M1=M+1 

M2=M+2 

M3=M+3 

go  to  (50,60,70) ,ki 
50  Z(L2)=Z(L2)+ZL(M1) 

Z(L3)=Z(L3)+ZL(M2) 

GO  TO  80 
60  Z(L1)=Z(L1)+ZL(M1) 

Z(L2)=Z(L2)+ZL(M2) 

Z(L3)=Z(L3)+ZL(M3) 

GO  TO  80 
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70  Z(L1)=Z(L1)+ZL(M1) 
Z(L2)=Z(L2)+ZL(M2) 
80  Z(L0+I)=Z(L0+I)+ZL(M0+I) 
100  CONTINUE 
RETURN 
END 
C*** ******************************** *************SUBROUTINE  OGIVE 


OGIVE 

4  SEPTEMBER  1991 
R.M.  FRANCIS 
26  JANUARY  1992 


C  SUBROUTINE 

C  DATE 

C  PROGRAMMER 

C  REVISED 

C  COMMENTS 

C  THIS  SUBROUTINE  WILL  GENERATE  DATA  FOR  A  BODY  OF  REVOLUTION  (BOR)  IN 

C  THE  FORM  OF  AN  OGIVE. 

C  DIMENSIONS  ARE  NORMALIZED  TO  WAVELENGTH,  SEE  PAGE  14  FRANCIS  THESIS. 

C  ZH  =  Z  CO-ORDINATE  *  2*PI 

C  RH  =  RADIUS  *2*PI 

SUBROUTINE  OGIVE ( NP , ARAD , ZH , RH , BASE , RS , ZP ) 
C  NP  =  NUMBER  OF  POINTS  ON  THE  OGIVE  SURFACE,  MAXIMUM  =  1000 
C  ZP   =  ZPRIME,  THE  POSITION  ON  Z  WHERE  THE  RADIUS  OF  CURVATURE  STARTS 
C  BASE  =  BASE  RADIUS 
C  RS   =  RADIUS  OF  CURVATURE  IN  THE  RZ , Z  PLANE  OF  THE  OGIVE 

INTEGER  I.NP.NPBASE 

REAL  RH(400) , ZH (400 ) ,ZP, BASE, RS .ZCOORD, RADIUS ,AL 

PI=3. 1415926 
C  INPUT  THE  VARIABLES  FOR  THE  OGIVE, ZP, B ,RS,NP 

WRITE(6,*)  'ENTER  SURFACE  CURVATURE  (wavelengths)' 

READ(5,*)  RS 

WRITE(6,*)  'ENTER  ZPRIME, WHERE  CURVATURE  STARTS  (wavelengths)' 

READ(5,*)  ZP 

WRITE(6,*)  'ENTER  BASE  RADIUS  (wavelengths)' 

READ(5,*)  BASE 
C  PERFORM  CALCULATIONS 

AL=SQRT(2 . *BASE*RS-BASE**2) 
ZMAX=  AL  +  ZP 
ANG=ASIN(AL/RS) 

L=ANINT((RS*ANG+ZP)*5) 
NP=2*L+1 

WRITE (6,*) 'NUMBER  OF  POINTS  IS:  \NP 
DZ=  ZMAX/FL0AT(NP-1) 
DO  10   1=1, NP 

ZCOORD=ZMAX-  FL0AT(I-1)*DZ 
IF(ZCOORD.Eq.O.)THEN 

ZCOORD=0. 0000000001 
ENDIF 
ZH(I)=2.*PI*ZC00RD 

IF  (ZCOORD. LE.ZP)  THEN 

RADIUS=BASE 
ELSE 

RADIUS=SqRT(RS**2-(ZC00RD-ZP)**2)+(BASE-RS) 
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IF (RADIUS. Eq.O. )THEN 

RADIUS=0. 0000000001 
ENDIF 
ENDIF 
RH(I)=2.*PI*RADIUS 
10  CONTINUE 
RETURN 
END 
O***************************************************SUBR0UTINE  GENEX 


C  SUBROUTINE 
C  DATE 
C  REVISED 
C  PROGRAMMER 


GENEX 

14  JANUARY  1992 
17  February  1992 
R.M.  FRANCIS 

SUBROUTINE  GENEX (MODE , NP , NT , NPHI , CNPHI , XT , AT , X , A , 
*  XP.AP.THS.PHS.ARAD.RHO.ZHO.R) 

C  SOURCE  ON  Z  AXIS  AT  Z=0.  THIS  VERSION  DOES  THE  ANTENNA  INTEGRATION 
C  IN  RECTANGULAR  COORDINATES.   CNPHI  POINTS  USED  IN  X  AND  Y . 
C  **************************** 

c  * 

C  t          t    i     * 

C  (V)     =<W,E>    * 

C  n   l       ni         * 

C  * 

C  * 

C  phi        phi   i     * 

C  (V)    =<W,E>    * 

C  n  i       ni        * 

C  * 
<3***************************** 

C  THE  EXCITATION  VECTOR  IS  COMPUTED  FOR  THE  GIVEN  R.THETA  AND  PHI 
C  COMPONENTS  ON  SURFACE  OF  BOR,  SPECIFIED  IN  SUBROUTINE  CIRCRTP. 

COMPLEX  R(800) ,CEXP,PSI,S1,S2,S3,S4,S5 

COMPLEX  CIRCR.CIRCT.CIRCP 

DIMENSION  RH(400) ,ZH(400) ,XT(4) ,AT(4) ,X(100) ,A(100) 

DIMENSION  XP(IOO) ,AP(100) ,RH0(400) ,ZH0(400) 

INTEGER  CNPHI 

PI=3. 141592654 

BK=2.*PI 

MP=NP-1 

MT=MP-1 

N=MT+MP 
C  LIMITS  ON  PHI  INTEGRATION 

Pl=(2.*PI-0.)/2. 

P2=(2.*PI+0.)/2. 

DO  10  1  =  1, NP 

RH(I)=RHO(I)/BK 
10    ZH(I)=ZHO(I)/BK 

DO  30  IP=1,MP 
C  QUANTITIES  FOR  THE  FIRST  SEGMENT  (POSITIVE  SLOPE) 

I=IP+1 

II=IP 


89 


DR=RH(I)-RH(II) 

DZ=ZH(I)-ZH(II) 

D1=SQRT(DR*DR+DZ*DZ) 

Rl=.5*(Rfl(I)+RH(II)) 

Z1=.5*(ZH(I)+ZH(II)) 

SVP1=DR/D1 

CVP1=DZ/D1 

Vl=ATAN2(DR,DZ+l.E-5) 
C  QUANTITIES  FOR  THE  SECOND  SEGMENT  (NEGATIVE  SLOPE) 
C  (SKIP  IF  LAST  SEGMENT) 

I=IP+2 

II=IP+1 

DR=RH(I)-RH(II) 

DZ=ZH(I)-ZH(II) 

D2=SQRT(DR*DR+DZ*DZ) 

R2=.5*(Rfl(I)+RH(II)) 

Z2=.5*(ZH(I)+ZH(II)) 

SVP2=DR/D2 

CVP2=DZ/D2 

V2=ATAN2(DR,DZ+1 .E-5) 
C  BEGIN  PHI  INTEGRATION:  R  TERMS  IN  SI  AND  S2;  THETA  TERM  IN  S3  AND  S4; 
C  PHI  TERM  IN  S5. 

S1=(0.,0.) 

S2=(0. ,0.) 

S3=(0. ,0.) 

S4=(0. ,0.) 

S5=(0. ,0.) 

DO  20  J=1,NPHI   - 

PH=P1*X(J)+P2 

PSI=CEXP(CMPLX(0. ,-MODE*PH))*A(J) 

IF(IP.LE.MT)  THEN 
C  t-CURRENT  INTEGRATION  FOR  THE  POSITIVE  SLOPE 
C  Gauss  Quadrature 

DO  13  L=1,NT 

TP=Dl*XT(L)/2. 

RHB=R1+TP*SVP1 

ZHB=Z1+TP*CVP1 

TH=ATAN2(RHB .ZHB+1 . E-5) 

CC=C0S(V1-TH) 

SS=SIN(V1-TH) 

CALL  CIRCRTP (CNPHI , XP , AP , ARAD , THS , PHS , 
*  PH.RHB.ZHB.CIRCR.CIRCT.CIRCP) 

S1=S1+AT(L)*(.5+TP/D1)*CC*CIRCR*PSI 

S2=S2+AT(L)*(.5+TP/D1)*SS*CIRCT*PSI 
13    CONTINUE 

C  t-CURRENT  INTEGRATION  FOR  THE  NEGATIVE  SLOPE 
C  Gauss  Quadrature 

DO  14  L=1,NT 

TP=D2*XT(L)/2. 

RHB=R2+TP*SVP2 

ZHB=Z2+TP*CVP2 
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TH=ATAN2(RHB , ZHB+1 . E-5) 

SS=SIN(V2-TH) 

CC=C0S(V2-TH) 

CALL  CIRCRTP ( CNPHI , XP , AP , ARAD , THS , PHS , 

*  PH.RHB.ZHB.CIRCR.CIRCT.CIRCP) 
S3=S3+AT(L)*(.5-TP/D2)*CC*CIRCR*PSI 
S4=S4+AT(L)*(.5-TP/D2)*SS*CIRCT*PSI 

14  CONTINUE 
ENDIF 

C  phi-CURRENT  INTEGRATION 
DO  15  L=1,NT 
TP=Dl*XT(L)/2. 
RHB=R1+TP*SVP1 
ZHB=Z1+TP*CVP1 
TH=ATAN2(RHB .ZHB+1 . E-5) 
CALL  CIRCRTP ( CNPHI , XP , AP , ARAD , THS , PHS , 

*  PH,RHB,ZHB,CIRCR,CIRCT,CIRCP) 
S5=S5+AT(L)*CIRCP*RHB/R1*PSI 

15  CONTINUE 
20    CONTINUE 

C  COMPONENTS  ARE  STORED  IN  A  COLUMN  VECTOR:  VT(1,MT),  VP(MT+1,N) 

IF(IP.LE.MT)  THEN 

R(IP)=(S1+S2)*D1/2.*P1+(S3+S4)*D2/2.*P1 

ENDIF 

R(IP+MT)=S5*Pl*Dl/2. 
30    CONTINUE 

RETURN 

END 
r***** *************** ******************************SUBROUTINE  CIRCRTP 


CIRCRTP 

1  OCTOBER  1991 

26  February  1992 
R.M.  FRANCIS 

27  October  1992 
K.  A.  KLOPP 


C  SUBROUTINE 

C  DATE 

C  REVISED 

C  PROGRAMMER 

C  2ND  REVISION 

C  PROGRAMMER 

C 

C  COMMENTS  :   THIS  SUBROUTINE  WILL  RIGOROUSLY  CALCULATE  THE  ELECTRIC  FIELD 

C  FOR  A  CIRCULAR  APERTURE.  THE  FIELD  IS  CALCULATED  AT  THE  COORDINATES 

C  SPECIFIED  BY  RH(.)  AND  ZH( . )   THE  APERTURE  IS  LOCATED  AT  Z  =  0 , 

C  AND  SCANNED  TO  A  DIRECTION  (THS, PHS).  THE  SUBROUTINE  OGIVE  IS  THE 

C  SOURCE  OF  THE  GEOMETRIC  DATA  REQUIRED  BY  CIRC  TO  PERFORM  COMPUTATIONS. 

C  ALL  PHYSICAL  DIMENSIONS  ARE  NORMALIZED  TO  WAVELENGTH. 

C 

c***************************************************************************** 

C                   +ant  +ant  * 

C  * 

C    l               \  \  xl  * 

C  E  =  -jn/(4*pi*k)   /_  /_   e  ( )*TAPER(rp)  * 

C                    x=0  y=0  * 

C  * 

I;****************************************************************************" 


SUBROUTINE  CIRCRTP (CNPHI , XP , AP , ARAD , THS , PHS , 
*  PHI,RZ,Z(CIRCR,CIRCT,CIRCP) 

INTEGER  CNPHI 

DIMENSION  XP(100),AP(100) 

REAL  K 

COMPLEX  J , JK  ,  Gl , G2 , XI , CON , CC , SUMX , SUMY , SUMZ 

COMPLEX  CIRCR.CIRCT.CIRCP 

PI=3. 141592654 

ETA=120.*PI 
C  SET  FIELD  COMPONENTS  TO  ZERO  IN  THE  REAR  HEMISPHERE 

CIRCR=(0. ,0. ) 

CIRCT=(0. ,0. ) 

CIRCP=(0. ,0. ) 
C  Z  OF  THE  OBSERVATION  POINT  IS  >=  0  IN  THE  FORWARD  HEMISPHERE 

IF(Z.GE.O.)  THEN 

SUMX=(0. 0,0.0) 

SUMY=(0. 0,0.0) 

SUMZ=(0. 0,0.0) 

K=6. 283185307 

J  =  (0.0,1.0) 

JK=  (0.0,6.283185307) 

C0N=-J/(4.*PI) 
C  OMIT  THE  FACTOR  OF  ETA  SINCE  IT  IS  NOT  INCLUDED  IN  SUBROUTINE  ZMAT. 
C  ARAD**2  IS  FOR  SCALING  OF  GAUSSIAN  INTEGRATION.   READJUST  SCALING 
C  TO  AGREE  WITH  OLD  SUBROUTINE  SPHERE  (ARBITRARY) 

CC=  CON 
C  OUTER  LOOP:  INTEGRATE  OVER  X,  -ARAD  <  X  <  ARAD. 
C  INNER  LOOP:  INTEGRATE  OVER  Y,  -ARAD  <  Y  <  ARAD. 
C  EVALUATE  ANTENNA  FIELD  AT  RZ.Z: 

XO=RZ*COS(PHI) 

YO=RZ*SIN(PHI) 
STS=SIN(-THS) 
S=RZ/SQRT(RZ**2+Z**2) 
C=Z/SqRT(RZ**2+Z**2) 
C  INTEGRATE  IN  X,Y  COORDINATES 
DO  50  M=l, CNPHI 
X=ARAD*XP(M) 
DO  60  N=l, CNPHI 
Y=ARAD*XP(N) 
RP=SQRT(X**2+Y**2) 
C  SKIP  INTEGRATION  POINTS  BEYOND  THE  ANTENNA  RADIUS 

IF(RP.LE.ARAD)  THEN 

PHIPRIME=ATAN2(Y ,X+1 . E-10) 

R=SQRT(((RZ-RP)**2+Z**2)+(4*RZ*RP*SIN((PHI-PHIPRIME)/2)**2)) 

G1=(((K*R)**2)-1-(JK*R))/R**3 

G2=(3+(3*JK*R)-(K*R)**2)/R**5 

X1=JK*(RP*C0S(PHIPRIME-PHS)*STS-R) 

Y1=X0-X 

X2=Y1**2 

Y2=Y0-Y 

SUMX=SUMX+(CC*AP(M)*AP(N)*(G1+(X2*G2))*CEXP(X1))*TAPER(RP) 
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SUMY=SUMY+(CC*AP(M)*AP(N)*Y1*Y2*G2*CEXP(X1))*TAPER(RP) 
SUMZ=SUMZ+(CC*AP(M)*AP(N)*Y1*Z*G2*CEXP(X1))*TAPER(RP) 
ENDIF 
60    CONTINUE 
50    CONTINUE 
C  ASSUME  AN  ELEMENT  FACTOR,  ELMT 
c     ELMT=SQRT(ABS(C)) 
ELMT=1 . 

CIRCT=(C*COS(PHI)*SUMX+C*SIN(PHI)*SUMY-S*SUMZ)*ELMT 
CIRCR=(S*COS(PHI)*SUMX+S*SIN(PHI)*SUMY+C*SUMZ)*ELMT 
CIRCP= ( -SIN (PHI ) *SUMX+COS (PHI ) *SUMY) *ELMT 
ENDIF 
RETURN 
END 
C**********************************************SUB ROUTINE  TESTSPHERE 


C  SUBROUTINE 
C  DATE 

C  PROGRAMMER 
C  COMMENTS 


TESTSPHERE 
21  FEBRUARY  1992 
R.  M.  FRANCIS 

GENERATES  A  SPHERE,  WITH  PLOTTED  POINTS 
C  AT  EQUAL  INTERVALS  IN  THETA . 

SUBROUTINE  TESTSPHERE ( NP , ZH , RH , BASE , RS , ZP ) 
INTEGER  HP, I 

REAL  ZH (400 ) , RH (400 ) , BASE , RS , ZP , SPRAD 
PI=3. 1415926 
BK=2.*PI 
ZP=0.0 

WRITE (6,*) 'ENTER  NUMBER  OF  POINTS  (HP):' 
READ(5,*)NP 

WRITE (6,*) 'ENTER  SPHERE  RADIUS' 
READ (5,*) SPRAD 
BASE=SPRAD 
RS=SPRAD 
DO  1241  1=1, HP 

ANGLE=PI*FL0AT(I-1)/(2.*FL0AT(NP-1)) 
ZH(I)=BK*SPRAD*COS (ANGLE) 
RH(I)=BK*SPRAD*SIN( ANGLE) 
IF(RH(I).EQ.O.)THEN 

RH(I)=0. 0000000001 
ENDIF 
1241   CONTINUE 
RETURN 
END 
C*** *************************************** ********SUBROUTINE  CONE 
C  GENERATES  THE  BOR  GEOMETRY  FOR  A  CONE. 

SUBROUTINE  CONE ( NP , ZH , RH , BASE , RS , ZP ) 

REAL  ZH(400) ,RH(400) ,BASE,RS,ZP 

REAL  HA,ANG,ZMAX,DZH 

INTEGER  NP 

PI=3. 1415926 

BK=2.*PI 

WRITE(6,*) 'ENTER  THE  CONE  HALF  ANGLE  IN  DEGREES' 
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READ (5,*) HA 

IF(HA.GT.90. )THEN 

WRITE(6,*) 'RANGE  OF  ANGLE:  0<  ANGLEOO ' 

WRITE ( 6 , * ) ' RANGE  EXCEEDED ' 

GOTO  10 

ENDIF 
ANG=HA*PI/180. 

WRITE(6,*) 'ENTER  THE  MAXIMUM  Z' 
READ(5,*)ZMAX 

WRITE (6,*) 'ENTER  THE  NUMBER  OF  POINTS  (NP) ' 
READ(5,*)NP 
ZH(1)=ZMAX*BK 
RH(1)=0. 

DZH=ZH(1)/FL0AT(NP-1.  ) 
DO  10  1=1, NP 

IF(ZH(I).EQ.O.)THEN 
ZH(I)=. 00000001 

ENDIF 
ZH(I+1)=ZH(I)-DZH 
RH(I+1)=RH(I)+(DZH*SIN(ANG)) 
10  CONTINUE 

BASE=RH(NP)/BK 
ZP=0. 
RS=0. 
RETURN 
END 
C* ************************************** *******SUBROUTINE  DISK 
SUBROUTINE  DISK (NP , ZH , RH , BASE , RS , ZP ) 
REAL  ZH(400) ,RH(400) ,BASE,RS,ZP ,RMAX 
INTEGER  NP 

WRITE(6,*) 'ENTER  THE  RADIUS  OF  THE  DISK' 
READ(5,*)RMAX 

WRITE (6,*) 'ENTER  THE  DISTANCE  TO  THE  DISK' 
READ(5,*)ZP 

WRITE (6,*) 'ENTER  THE  NUMBER  OF  POINTS  (NP)' 
READ(5,*)NP 
PI=3. 1415926 
BK=2.*PI 

DRH=RMAX*BK/FL0AT(NP-1 . ) 
DO  10  1=1, NP 

ZH(I)=ZP*BK 

RH(I)=DRH*FL0AT(I-1.) 
10  CONTINUE 
BASE=RMAX 
RS=0. 
RETURN 
END 
C* **************************** *****************SUBR0UTINE  PARAB 
SUBROUTINE  PARAB(NP ,ZH,RH,DM,FOD ,ZP) 
DIMENSION  RH(400),ZH(400) 
WRITE (6,*)  'ENTER  DIAMETER  AND  F/D  RATIO' 
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READ(5,*)  DM.FOD 

WRITE(6,*)  'FEED  DISTANCE  FROM  FOCUS  (+  IS  NEARER)' 

READ(5,*)  ZP 

WRITE(6,*)  'ENTER  NUMBER  OF  GENERATING  POINTS' 

READ(5,*)  NP 

BK=2.*3. 14159 

FM=FOD*DM 

PHIV=2.*ATAN(l./4./F0D) 

DO  24  1=1, HP 

TH=FLOAT(I-l)*PHIV/FLOAT(NP-l) 

RM=2. *FM/(l .+COS(TH) )*BK 

ZH(I)=RM*COS(TH)+ZP*BK 

RH(I)=RM*SIN(TH) 
24    CONTINUE 

RETURN 

END 
C********** *********************** *********FUNCTION  TAPER 

FUNCTION  TAPER(RHO) 
C  SPECIFY  ANTENNA  ILLUMINATION  FUNCTION.   REAL  FUNCTION  OF 
C  RHO  ONLY 

TAPER=1. 

RETURN 

END 
C* ***************** ***********************SUBR0UTINE  ANTFF 

SUBROUTINE  ANTFF(THETA , PHI , THETAS , PHIS , ARAD , ETF , EPF) 
C  COMPUTES  CLOSED  FORM  ANTENNA  PATTERN  IN  THE  FAR  FIELD. 
C  ETF  AND  EPF  ARE  RETURNED.  ANGLES  ARE  IN  RADIANS. 
C  PRESENT  ENTRY  IS  FOR: 
C 

C  UNIFORM  ILLUMINATION  SCANNED  TO  A  DIRECTION  (THETAS, PHIS) 
C 

COMPLEX  ETF,EPF,JK,SCL,EB 

JK=(0. 0,6. 283185307) 

BK=6. 283185307 

PI=BK/2. 
C  SET  RADIATION  PATTERN  TO  ZERO  IN  THE  REAR  HEMISPHERE 

IF(COS(THETA).LT.O.)  THEN 
EB=(0. ,0.) 

ELSE 

BB=SIN(THETA)*SIN(PHI)-SIN(THETAS)*SIN(PHIS) 
AA=SIN(THETA)*COS(PHI)-SIN(THETAS)*COS(PHIS) 
ARG=BK*ARAD*SQRT(AA**2+BB**2) 
SCL=-(0. ,1.)*2.*PI**2 
IF(ABS(ARG) .LT.1.E-5)THEN 
EB=SCL/2. 
ELSE 

EB=SCL*(BESSJ1(ARG)/ARG) 
ENDIF 

ENDIF 
C  EB  IS  THE  X  COMPONENT;  GET  ETHETA  (ETF)  AND  EPHI  (EPF) 
C  ASSUME  AN  ELEMENT  FACTOR,  ELMT 
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ELMT=SqRT(ABS(COS(THETA) ) ) 
ELMT=1. 
ETF=EB*COS (THETA ) *COS (PHI ) *ELMT 
EPF=-EB*SIN (PHI ) *ELMT 
RETURN 
END 

FUNCTION  BESSJ(N.X) 
C  RETURNS  THE  BESSEL  FUNCTION  B  OF  ORDER  N  (>1)  AND  REAL 
C  ARGUMENT  X. 

PARAMETER  (IACC=40,BIGNO=l . ElO ,BIGNI=1 . E-10) 
IF(N.EQ.O)  THEN 
C  IF  N=0  CALL  BESSJO 

BESSJ=BESSJO(X) 
ELSE  IF(N.EQ.l)  THEN 
C  IF  N=l  CALL  BESSJ1 

BESSJ=BESSJ1(X) 
ELSE 
C  IF  N>1  USE  RECURSION 
BESSJ=0. 

IF(X.NE.O.)  THEN 
T0X=2./X 

IF(X.GT.FLOAT(N))  THEN 
BJM=BESSJO(X) 
BJ=BESSJ1(X) 
DO  11  J=1(N-1 
BJP=J*TOX*BJ-BJM 
BJM=BJ 
BJ=BJP 

11  CONTINUE 
BESSJ=BJ 
ELSE 

M=2*( (N+INT(SQRT(FLOAT(IACC*N) ) ) )/2) 

BESSJ=0. 

JSUM=0. 

SUM=0. 

BJP=0. 

33-1. 

DO  12  J=M,1,-1 

BJM=J*TOX*BJ-BJP 

BJP=BJ 

BJ=BJM 

IF(ABS(BJ).GT.BIGNO)  THEN 

BJ=BJ*BIGNI 

BJP=BJP*BIGNI 

BESSJ=BESSJ*BIGNI 

SUM=SUM*BIGNI 

ENDIF 

IF(JSUM.NE.O)  SUM=SUM+BJ 

JSUM=1-JSUM 

IF(J.EQ.N)  BESSJ=BJP 

12  CONTINUE 


96 


SUM=2.*SUM-BJ 

BESSJ=BESSJ/SUM 

ENDIF 

ENDIF 

ENDIF 

RETURN 

END 

FUNCTION  BESSJO(X) 
C 

C  BESSEL  FUNCTION  OF  0  ORDER,  REAL  ARGUMENT  X 
C  (SEE  'NUMERICAL  RECIPES',  P. 172) 
C 

REAL*8  Y , P 1 , P2 , P3 , P4 , P5 , q 1 , q2 , Q3 , Q4 , q5 , Rl , R2 , R3 , R4 , R5 , R6 , 

*  S1,S2,S3,S4,S5,S6 

DATA  P1,P2,P3,P4,P5/1.D0,-. 109862827D-2 , . 2734510407D-4, 

*  - . 2073370639D-5 , . 209388721 1D-6/ 

DATA  Ql,Q2,q3,Q4,q5/-. 1562499995D-1, . 1430488765D-3 , 

*  -.6911147651D-5, . 7621095161D-6 ,- . 934945152D-7/ 

DATA  Rl , R2 , R3 , R4 , R5 , R6/57568490574 . DO , - 13362590354 . DO , 

*  651619640 . 7D0 , - 1 1214424 . 18D0 , 77392 . 33017D0 , -184 . 9052456D0/ 
DATA  S 1 , S2 , S3 , S4 , S5 , S6/575684904 1 1 . DO , 1029532985 . DO , 

*  9494680 . 718D0 , 59272 . 64853D0 , 267 . 8532712D0 , 1 . DO/ 
IF(X.EQ.O.)  THEN 

BESSJ0=1. 
ELSE  IF(ABS(X).LT.8.)  THEN 
Y=X**2 
BESSJ0=(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6)))))/ 

*  (S1+Y*(S2+Y*(.S3+Y*(S4+Y*(S5+Y*S6))))) 
ELSE 

AX=ABS(X) 

Z=8./AX 

Y=Z**2 

XX=AX-. 785398164 

BESSJ0=SqRT(.636619772/AX)*(C0S(XX)*(Pl+Y*(P2+Y*(P3+ 

*  Y*(P4+Y*P5))))-Z*SIN(XX)*(qi+Y*(q2+Y*(q3+ 

*  Y*(q4+Y*q5))))) 

ENDIF 

RETURN 

END 

FUNCTION  BESSJl(X) 
C 

C  BESSEL  FUNCTION  B  OF  ORDER  1,  REAL  ARGUMENT  X 
C  (SEE  'NUMERICAL  RECIPES',  P. 173) 
C 

REAL*8  Y , PI , P2 , P3 , P4 , P5 , qi , q2 , q3 , q4 , q5 , Rl , R2 , R3 , R4 , R5 , R6 , 

*  S1,S2,S3,S4,S5,S6 

DATA  P1,P2,P3,P4,P5/1.D0, . 183105D-2 , - . 3516396496D-4 , 

*  . 2457520174D-5.-.20337019D-6/ 

DATA  qi , q2 , q3 , q4 , q5/ . 04687499995D0 , - . 2002690873D-3 , 

*  . 8449199096D-5 , - . 99228987D-6 , . 105787412D-6/ 

DATA  Rl , R2 , R3 , R4 , R5 , R6/723626 14232 . DO , -7895059235 . DO , 
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*  242396853 . 1D0 , -297261 1 . 439D0 , 15704 . 4826D0 , -30 . 16036606D0/ 
DATA  SI , S2 , S3 , S4 , S5 , S6/144725228442 . DO , 2300535178 . DO , 

*  18583304 . 74D0 , 99447 . 43394D0 , 376 . 9991397D0 , 1 . DO/ 
IFCX.EQ.O. )  THEN 

BESSJ1=0. 
ELSE  IF(ABS(X).LT.8.)  THEN 
Y=X**2 
BESSJ1=X*(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6)))))/ 

*  (S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6))))) 
ELSE 

AX=ABS(X) 

Z=8./AX 

Y=Z**2 

XX=AX-2. 356194491 

BESSJ1=SQRT(.636619772/AX)*(C0S(XX)*(P1+Y*(P2+Y*(P3+ 

*  Y*(P4+Y*P5) ) ) )-Z*SIN(XX)*(qi+Y*(q2+Y*(Q3+ 

*  Y*(Q4+Y*Q5)))))*SIGN(1. ,X) 
ENDIF 

RETURN 
END 
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APPENDIX  C.  GAIN.F  LISTING 

Material  on  the  following  pages  is  a  listing  of  the  program  GAIN.F  described 
in  chapter  III. 
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c*********************************************************************** 
c 

C  PROGRAM     :  gain.f 

C  DATE        :  17  November  1992 

C  PROGRAMMERS  :  Dr.  D.  Jenn,  LT  K.  Klopp 

C 

C  Computes  the  gain  of  a  circular  aperature  antenna  radiating  through 

C  a  radome  consisting  of  a  body  of  revolution,  by  numerically 

C  integrating  electric  field  surrounding  system.   E-Field  is  generated  by 

C  currents  distributed  on  a  surface  of  the  body  of  revolution.   The 

C  currents,  and  their  locations  were  previously  determined  using  the 

C  program  "radome.f",  and  stored  in  the  file  "f77curcoef sA  or  B." 

C 

C  INPUT  <  f77curcoefsA/B 

C 

C*********************************************************************** 

CHARACTER  ITH 

CHARACTER* 12  CC 

CHARACTER* 10  GO 

CHARACTER*8  GNT.GNPHI ,CGP 

CHARACTER*14  TPTS ,PPTS .PHIPTS 

COMPLEX  R(1600) , LEADTERM , ETT , EPT 

COMPLEX  EP , ET , ETF , EPF , ETCOMB , EPCOMB 

COMPLEX  EXP1,EXP2,C0NJG,CEXP,CMPLX,CT(30000) .IMP.JK 

COMPLEX  U,UC,ET1,EP1,ET2,EP2,ZL0(400) 
C      DIMENSION  RH(400),ZH(400),IPS(800) 

DIMENSION  RH(400),ZH(400) 

DIMENSION  XT(4) ,AT(4) , A(100) ,X( 100) ,XP( 100) ,AP(100) 

DIMENSION  Q(30),S(30) 

DIMENSION  THRR(IOOO.IOOO) ,PHRR(1000) 
C      DIMENSION  ETSCAT(500),EPSCAT(500) 
C      DIMENSION  EXP(500),ANG(500),ECP(500) 

INTEGER  CNPHI, SELECTION 

DATA  PI/3.1415926/ 

DATA  RR/1000/ 

Rad=PI/180. 

BK=2.*PI 

U=(0. ,1.) 

U0=(0. ,0. ) 

UC=-U/4./PI 

JK=(0. 0,6. 283185307) 
C*****ENTER  A  LETTER  TO  INDICATE  WHICH  ITERATION  THIS  IS***************** 
C*****    (ALL  DATA  FILES  ARE  APPENDED  WITH  THIS  LETTER)    *********** 
C*****    (.m  FILES  ARE  FOR  GENERATING  PLOTS  IN  MATLAB)    *********** 

WRITE (6,*) 'ENTER  A  LETTER  TO  INDICATE  WHICH  ITERATION  THIS  IS' 

READ(5,*)ITH 

CC= ' curcoef sdat ' //ITH 

G0= ' gainout ' //ITH// ' . m ' 
C************************************************************************ 

c 

C   READ  in  the  following  data  from  file  called  f77curcoefs: 
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c 

C         iteration,  ITH 

C        geometry  of  BOR,  SELECTION 

C        radius  of  base  of  BOR,       BASE 

C        total  points,  NP 

C         number  of  t- intervals,       MT 

C        phi-intervals ,  MP 

C         sum  of  t  and  phi  intervals,   N 

C         number  of  rows  of  currents,   NROW 

C         total  number  of  modes ,       MODES 

C         number  of  sub  vectors,       NBLOCK 

C         antenna  scan  angle  in  theta,  THS 

C         antenna  scan  angle  in  phi,    PHS 

C        antenna  radius,  ARAD 

C        observation  angle,  PHIO 

C         complex  radome  impedance,     IMP 

C        zprime,  ZP 

C        surface  radius,  RS 

C         closed  form  antenna  pattern,  IFF=0 

C 

C        files  where  Gauss  integration  weights  and  abscissas  are  stored 

C 

C        coordinates  of  surface  points  *2PI/lamda,  ZH(I),  and  RH(I) 

C 

C ************************************************************************* 

OPEN(l,file=CC) 

WRITE(6,*)  'FILE  OPENED  IS  \CC 

READ ( 1 , * )  ITH , SELECTION , BASE , NP , MT , MP , N 

READ ( 1 , * )  NROW , MODES , NBLOCK , THS , PHS , ARAD , PHIO 

READ(1,*)  IMP.ZP.RS.IFF 
c      write(6,*)IMP,ZP,RS,IFF 

READU,*)  (RH(L),L=1,NP) 
c      write(6,*)  (RH(L) ,L=1 ,NP) 

READ(1,*)  (ZH(L),L=1,NP) 
c      write(6,*)  (ZH(L) ,L=1 ,NP) 

READ(1,*)  (ZL0(L),L=1,NP) 
C      write(6,*)  (ZLO(L) ,L=1,NP) 
C ************************************************************************** 

c 

C   READ  in  current  coefficients,  CT(I)  from  file  curcoef sdatA/B 

C    (This  is  one  one  collumn  vector.) 

C 

READ(1,*)  (CT(L),L=1,NR0W) 
C       wnte(6,*)(CT(L),L=l,NR0W) 
READ(1,*)  GNT.GNPHI.CGP 

WRITE(6,*)  'RUN  CODE  READ  FROM  DISC  :  ',ITH 
THETAS=THS*RAD 
PHIS=PHS*RAD 
MHI=M0DES+1 
CLOSE(l) 


c 

C   Read  in  weights  A(I),  and  Abscissas  X(I)  for  the  following  Gauss- 

C   quadrature  Integrations: 

C  file      wts  absc 

C 

C      BOR  integration  in  t  direction:        TPTS   =  gaus/gaus#;  XT(I),AT(I) 

C      BOR  integration  in  phi  direction:      NPHI   =  gaus/gaus#;  X  (I), A  (I) 

C      Antenna  integration  in  x,y  direction:   PHIPTS  =  gaus/gaus#;  XP(I),AP(I) 

C 

TPTS='gaus/'//GNT 
PPTS='gaus/V/GNPHI 
PHIPTS='gaus/'//CGP 
OPEN ( 2 , FILE=TPTS , STATUS= ' OLD ' ) 
0PEN(3,FILE=PPTS,STATUS='0LD') 
OPEN ( 8 , FILE=PHIPTS , STATUS= ' OLD ' ) 
READ(2,*)NT 
WRITE(6,«) 'NT  =  ' ,NT 
IF(NT.GT.4)THEN 

WRITE (6,*) 'MAXIMUM  NUMBER  OF  POINTS(NT)  IS  4' 
GOTO  999 
ENDIF 
READ(3,*)NPHI 
WRITE(6,*)  'NPHI  =  \NPHT 
IF(NPHI.GT.200)THEN 

WRITE(6,*) 'MAXIMUM  NUMBER  OF  POINTS(NPHI)  IS  200' 
GOTO  999 
ENDIF 
READ (8,*) CNPHI 
WRITE (6,*) 'CNPHI  =  \CNPHI 
IF(CNPHI.GT.100)THEN 

WRITE (6,*) 'MAXIMUM  NUMBER  OF  POINTS(CNPBI)  IS  100' 
GOTO  999 
ENDIF 
C  LOAD  THE  WEIGHTS  AND  ABSCISSAS  IN  THE  VECTORS. 
DO  1  M=1,NT 

READ(2,*,END=1)XT(M) ,AT(M) 
WRITE(6,*)XT(M),AT(M) 

1  CONTINUE 

DO  2  M=1,NPHI 

READ(3,*,END=2)X(M),A(M) 
WRITE(6,*)X(M),A(M) 

2  CONTINUE 

DO  4  M=l, CNPHI 

READ(8,*,END=4)XP(M),AP(M) 

WRITE(6,*)XP(M),AP(M) 

4   CONTINUE 

CL0SE(2) 

CL0SE(3) 

CL0SEC8) 
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ft, 
ft, 

WRITE(7,2001) 
ft, 


CI****************************************************************************** 

c 

C   Create  an  output  file  called  'gainoutA/B .m' ,and  output  data  on  geometry  etc. 

C 

,  7  +  ****, *«««**.««**.«**********♦*****»*«  +  *♦*  +  «**»**»**«»***.**««»***»*. ...,.».,. 

0PEN(7,FILE=G0) 

WRITE(7,2000)'  GEOMETRY  =  ' .SELECTION , '  NP  =  ' , NP 

MT  =  * ,MT, '  MP  =  '  ,MP,  '  N  =  ',H 

MODES  =  ' .MODES, '  NBLOCK  =  '.NBLOCK 

PHI  =  ' ,PHI 

BASE  =  '.BASE,'  ARAD  =  ' , ARAD , '  THETASCAN  =  ' ,THS 

PHISCAN  =  * ,PHS 
WRITE(7,*) 'NPHI  =' ,NPHI 

2000  FORMAT ( //5X,  ' y//.  RADOME  GAIN  CALCULATION  '/,'/,'// ,10 (/ "/,    '  ,A13,I4)) 

2001  FORMAT  ( 3  (/"/.  '  , A13 ,F9 . 2) ,/) 

C*************** ******************************************** ******************** 

c 

C   Gain  Computation: 

C  4  pi  *  Max  Radiation  Intensity 

C  G  = 

C  Integration  Over  Closed  Surface  (  Radiation 

C  Intensity  *  Beam  Solid  Angle  ) 

C 

C   Radiation  Intensity  =  (l/(2*eta) )* I  E(t,phi,r)  1**2 

C 

C   Since  E  is  in  the  far  field,  need  to  only  determine  Etheta,  Ephi . 

C 

C   |E|**2  =  IE  scatteredl**2  +  IE  antenna|**2 

C 

C   IE  scattered  I  =  Etheta  +  Ephi 

C 

Q******************************************************************************* 

C**Enter  #  divisions  for  integration  in  theta,  and  phi 

C  directions.  (This  allows  easy  variation  of  the  number  of  integration  steps 
C  without  having  to  recompute  a  different  set  of  weights  and  abscissas  for 
C   Gauss  Integration.) 

WRITE (6,*) 'ENTER  NUMBER  OF  DIVISIONS  IN  THETA  DIRECTION', 
&  '  NDIVT  =  ?' 

READ(5,*)NDIVT 

WRITE(6,*) 'NDIVT  =  ', NDIVT 

WRITE (7,*)"/.  NDIVT  =  ', NDIVT 

WRITE (6,*) 'ENTER  NUMBER  OF  DIVISIONS  IN  PHI  DIRECTION,  NDIVP  =  ?' 

READ(5,*)NDIVP 

WRITE(6,*) 'NDIVP  =  ', NDIVP 

WRITE(7, *)"/.  NDIVP  =  ', NDIVP 

kk=NDIVT 
C      do  500  HDIVT=l,kk 

NDIVP=NDIVT 
C***********************************integration  intervals  0<theta<180,  0<phi<360 

THETA 1=0 

THETA2=180 


103 


S1=THETA1 
S2=THETA2 
Sl=Sl*Rad 
S2=S2*Rad 
PHI1=0 
PHI2=360 
Q1=PHI1 
Q2=PHI2 
Ql=qi*Rad 
q2=Q2*Rad 
C***********************************theta  integration  angles 
DS=(S2-S1)/FL0AT(NDIVT) 
DO  200  I=1,NDIVT+1 
S(I)=(I-1)*DS 

200  CONTINUE 
C***********************************phi  integration  angles 

DQ=(Q2-qi)/FL0AT(NDIVP) 
DO  201  I=1,NDIVP+1 

q(i)=(i-i)*Dq 

201  CONTINUE 
SUM=0. 
SUMN0NE=0. 
UMAX=0. 
UMAXN0NE=0. 
EMAX=0. 
EMAXN0NE=0. 

C************************V  phi  integration  of  rad  intensity  *  beam  solid  angle  V 
DO  300  JJ=1,NDIVP 
Pl=Dq/2. 

P2=(q(jj+D+q(jj))/2. 

C**********************************v  phi  integration  determining  Escattered  V 
C  V   within  each  interval  of  phi,  NDIVP   V 

DO  300  J=1,NPHI 
PHR=P1*X(J)+P2 
C*****************V  theta  integration  of  rad  intensity  *  beam  solid  angle  V 
DO  301  II=1,NDIVT 
Tl=DS/2. 

T2=(S(II+l)+S(II))/2. 
C**************************V  theta  integration  determining  Escattered  V 
C  V   within  each  interval  of  theta,  NDIVT   V 

DO  301  I=1,NPHI 
THR=T1*X(I)+T2 
IF(THR.LT.PI)  GO  TO  302 
THR=(BK-THR) 
PHR=PHR+PI 
302         CONTINUE 
C***************************************************Beam  Solid  Angle 
STRA=SIN(THR)*A(I)*A(J) 
KJ=(JJ-1)*NPHI+J 
KI=(II-1)*NPHI+I 
THRR(KJ,KI)=THR*180/PI 
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PHRR(KJ)=PHR*180/PI 

SNP=SIN(PHR) 

CNP=COS(PHR) 

C **»*»*«***«»******♦***«******♦****«**«*****  +  »«+******♦****♦♦»**,»*»*,  ********* 

C 

c 

_  t.theta  _  phi.theta 

R  R 

n  n 


-jkr     MHI 


jnphi 


C 
C 
C        I      t      I  -jhe 

CI  I    =      \ 

C        I        si  4pir  / I    _  t,phi  _  phi, phi  I    I      _  phi 

C        I    E        I  n=l    I    R  R  III 

C        I      phi  I  n  n  n 

C 

CV********************************see  page  25.  H,E,  and  Combined  Field********V 

C  -jkr 

C  -jhe 

C  Leading  Term     :  eta  (h)  is  cancelled  by  its  inverse  in  CT  (I), 

C  4pir      r  dependence  is  ignored,  and  the  whole  Term  is 

C  scaled  by  2pi  for  gauss  quadrature  integration 

C 

C  Leading  Term  is  reduced  to  (-j/4pi)/2pi 

r ******»****♦*********************»***************,.**********»**************** v 

C  LEADTERM=UC/BK 

LEADTERM=UC 

ET=(0,0) 

EP=(0,0) 

ET1=(0.,0.) 

EP1=(0.  ,0.  ) 

ET2=(0. ,0.) 

EP2=(0.  ,0.) 
£**********»*♦*******»************************************»***  **************** y 

C  MHI 

C 

C  \ 

C  /__ 

C  n=l 

DO  305  M=1,MHI 
NM=M-1 

C*************************************************************  jnphi 
C  e 

EXP1=CEXP(CMPLX(0. ,FL0AT(NM)*PHR) ) 

EXP2=C0NJG(EXP1) 
C*************determme  plane  save  excitation  vector  in  far  field,  R 

CALL  PLAHE (NM , NM , NP , NT , RH , ZH , XT , AT , THR , R) 

NT0P1=M0DES-NM 

NT0P2=NBL0CK-(NT0P1+1) 

NS2=NT0P1*N 

NS1=NT0P2*N 
£*♦*****«*♦********************************************************* 
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C  I    _   t.thetal     I    _   t    I      jnphi  I    _   t.phi      I     I    _   t    I      jnphi 

C  IR  l.llle  , and      I    R  I .  I    I        I    e 

C  In  I    |      n      I  In  llnl 

C******************************************************************* 

DO  306  L=1,HT 
C*** ******************************************* *pos it ive  mode  +n 

ET1=ET1+R(L)*CT(L+NS1)*EXP1 

EP1=EP1+R(L+N)*CT(L+NS1)*EXP1 
C* *************************************** *******negat ive  mode  -n 

IF(NM.Eq.O)  GO  TO  306 

ET1=ET1+R(L)*CT(L+NS2)*EXP2 

EP1=EP1-R(L+N)*CT(L+NS2)*EXP2 

306  CONTINUE 

C******************************************************************* 
C      I    _  phi.thetal    I    _   t    I      jnphi  I    _  phi, phi      I    I    _  t    I      jnphi 

CIR  l.llle  ,andlR  l.llle 

Cln  llnl  In  llnl 

C******************************************************************* 

DO  307  L=1,MP 

C***********************************************pOS1tive  mode   +n 
ET2=ET2+R(L+MT)*CT(L+NS1+MT)*EXP1 
EP2=EP2+R(L+MT+N)*CT(L+NS1+MT)*EXP1 

C*** ***************** ***************************negat ive  mode  -n 
IF(NM.EQ.O)  GO  TO  307 
ET2=ET2-R(L+MT)*CT(L+NS2+MT)*EXP2 
EP2=EP2+R(L+MT+N)*CT(L+NS2+MT)*EXP2 

307  CONTINUE 

c  ET=ET+LEADTERM*(ET1+ET2) 

c  EP=EP+LEADTERM*(EP1+EP2) 

305  CONTINUE 

C*********************************************** ~  end  mode  summation  " 

ETT=LEADTERH*(ET1+ET2) 

EPT=LEADTERM*(EP1+EP2) 
******************************************add  in  antenna  contribution 
C  FEED  CONTRIBUTION  IN  THE  FAR  FIELD  IS  ETF.EPF 

C****************************************************BESSJ1 _ 

C  USE  CLOSED  FORM  FEED  EXPRESSION  FROM  SUBROUTINE  ANTFF  IF 
C  IFF=0;  OTHERWISE  USE  BRUTE  FORCE  EVALUATION  FROM  CIRCRTP 
C 

R0BS=1000. 
IF(IFF.EQ.O)  THEN 
CALL  ANTFF (THR.PHR.THS ,PHS , ARAD ,ETF,EPF) 
C       write(6,*)'  ET  =  ' .CABS(ET),'  ETF  =  '.CABS(ETF) 
ELSE 
RHB=ROBS*SIN(THR) 
ZHB=ROBS*COS(THR) 

CALL  CIRCRTP (CNPHI , XP , AP , ARAD , THS , PHS , 
&  PHR,RHB,ZHB,CIRCR,CIRCT(CIRCP) 

C  REMOVE  THE  1/R  DEPENDENCE  BECAUSE  EXP(-jkR)/R  IS  OMITTED  IN 
C  THE  SCATTERED  FIELDS  ET  AND  EP 
ETF=CIRCT*ROBS 
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EPF=CIRCP*ROBS 
C       wnte(6,*)'  ET  =  '  ,CABS(ET),'  ETF(ROBS)  =  \CABS(ETF) 

ENDIF 
C**********************************************total  E  theta  and  E  phi 
C       wnte(6,*)  'THR(\jj,j,ii,i,NM,  ')  =  '  ,  180*THR/pi 
ETCOMB=ETT+ETF 
EPCOMB=EPT+EPF 

EHAG=SQRT(CABS(ETC0MB)**2+CABS(EPC0MB)**2) 
EMAGN0NE=SQRT(CABS(ETF)**2+CABS(EPF)**2) 
IF(EMAG.GT.EMAX)THEN 
EMAX=EMAG 
EMAXNONE=EMAGNONE 

WRITE(6,*)*  EMAG(',J,I, ')  =  '  ,EMAG 
END  IF 
C****************************************determine  max  field  intensity 
UMAG^EMAG**2 
UMAGN0NE=EMAGN0NE**2 
IF (UMAG.GT. UMAX) THEN 
UMAX=UMAG 

THETAMAX= 180*THR/PI 
PHIMAX=180*PHR/PI 
C  WRITE(6,*)'  UMAX  =  '.UMAX 

C  WRITE (6,*)'  THETAMAX  =  \THETAMAX 

C  WRITE(6,*)'  PHIMAX  =  '.PHIMAX 

END  IF 
IF (UMAGNONE . GT . UMAXNONE ) THEN 

UMAXNONE=UMAGNONE 
END  IF 
C***********************muitiply  integration  by  beam  solid  angle  (STRA) 
SUM=SUM+UMAG*STRA 
SUMNONE=SUMNONE+UMAGNONE*STRA 
301      CONTINUE 

C**************************************************  end  theta  integration  " 
300   CONTINUE 

C********************************************************~  end  phi  integration  " 
C**************************scale  each  summation  by  (b-a)/2  for  Gauss  integration 
PRAD=T1*P1*SUM 
PRADN0NE=T1*P1*SUMN0NE 

WRITE(6,*) 'PRAD  =  \PRAD 
WRITE(6,*)'PRAD  NO  RADOME  =  '.PRADNONE 
C******************************************************************************* 
C   End  Integration  over  closed  surface  radiation  intensity  *  beam  solid  angle 

Qm ******************************************************************  compute   gain 
GAIN=4*PI*UMAX/PRAD 
GAINN0NE=4*PI*UMAXN0NE/PRADN0NE 
GDB=10.*aloglO(GAIN) 
GDBN0NE=10 . *aloglO(GAINNONE) 

WRITE(6,*) 'GAIN' .NDIVT, '  =  '.GAIN,'  IN  DB  =  ' ,GDB 
WRITE(6,*)'GAINN0NE  =  '.GAINNONE,'  IN  DB  =  '.GDBNONE 
WRITE(7,*) ' NDIVT (', NDIVT, ')=', NDIVT, ';  NDIVPC .NDIVP, ')=',NDIVP 
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WRITE(7,2003) 'GAINO ,  NDIVT,  ')='  .GAIN, ' ;GAINDB(' ,  NDIVT, ')='  ,GDB 

WRITE(7(2003)'GAINN0NE(' .NDIVT, ')  =  \GAINNONE, 
ft' ;GAINNONEDB(\NDIVT, ')   =  '  .GDBNONE 
2003   F0RMAT(/(A,I3,A,F12.5,A,I3,A,F10.5)) 
C500   continue 

STOP 
999  end 
C************************************************************END  OF  MAIN  PROGRAM 
C*** ****************************** *******************SUBR0UTINE  PLANE 
C  REFERENCE:  TECHNICAL  REPORT  TR-80-1  ; (pages  57-64) 

SUBROUTINE  PLANE ( Ml , M2 , NP , NT , RH , ZH , XT , AT , THR , R ) 
C 

C  PLANE  WAVE  EXCITATION  VECTOR  IN  THE  FAR  FIELD.  FROM  MAUTZ  AND  HARRINGTON. 
C  MODIFIED  TO  DO  ONLY  ONE  ANGLE  AND  FREQUENCY  PER  CALL. 
C  **  NEW  BESSEL  FUNCTION  ROUTINE  FROM  NUM.  RECIPES  ** 
C 

COMPLEX  R(1600) ,U,U1 .UA.UB ,FA( 1500) ,FB(1500) ,F2A ,F2B,F1A,F1B 

COMPLEX  U2,U3,U4,U5 

DIMENSION  RH(400) ,ZH(400) ,XT(4) , AT(4) ,R2(4) ,Z2(4) 

MP=NP-1 

MT=MP-1 

N=MT+MP 

N2=2*N 

CC=COS(THR) 

SS=SIN(THR) 

U=(0.,1.) 

U1=3.141593*U**M1 

M3=M1+1 

M4=M2+3 

IF(Ml.EQ.O)  M3=2 

M5=Ml+2 

M6=M2+2 

DO  12  IP=1,MP 

K2=IP 

I=IP+1 

DR=.5*(RH(I)-RH(IP)) 

DZ=.5*(ZH(I)-ZH(IP)) 

D1=SQRT(DR*DR+DZ*DZ) 

R1=.25*(RH(I)+RH(IP)) 

IFCR1.EQ.0.)  Rial. 

Z1=.5*(ZH(I)+ZH(IP)) 

DR=.5*DR 

D2=DR/R1 

DO  13  L=1,NT 

R2(L)=R1+DR*XT(L) 

Z2(L)=Z1+DZ*XT(L) 
13  CONTINUE 

D3=DR*CC 

D4=-DZ*SS 

D5=D1*CC 

DO  23  M=M3,M4 
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FA(M)=0. 

FB(M)=0. 

23  CONTINUE 

DO  15  L=1,NT 
X=SS*R2(L)*2. 
ARG=Z2(L)*CC 

UA=AT(L)*CMPLX(COS(ARG),SIN(ARG)) 
UB=XT(L)*UA 
C  THIS  LINE  REPLACES  HARRINGTON'S  BLOCK 
DO  25  M=M3,M4 
BES=BESSJ(M-2,X) 
FA(M)=BES*UA+FA(M) 
FB(M)=BES*UB+FB(M) 

25  CONTINUE 
15  CONTINUE 

IF(Ml.NE.O)  GO  TO  26 

FA(1)=-FA(3) 

FB(1)=-FB(3) 

26  UA=U1 

DO  27  M=M5,M6 

M7=M-1 

M8=M+1 

F2A=UA*(FA(M8)+FA(M7)) 

F2B=UA*(FB(M8)+FB(M7) ) 

UB=U*UA 

F1A=UB*(FA(M8)-FA(M7)) 

F1B=UB*(FB(M8)-FB(M7)) 

U4=D4*UA 

U2=D3*F1A+U4*FA(M) 

U3=D3*F1B+U4*FB(M) 

U4=DR*F2A 

U5=DR*F2B 

K1=K2-1 

K4=K1+N 

K5=K2+N 

R(K2+MT)=-D5*(F2A+D2*F2B) 

R(K5+MT)=D1*(F1A+D2*F1B) 

IF(IP.EQ.l)  GO  TO  21 

R(K1)=R(K1)+U2-U3 

R(K4)=R(K4)+U4-U5 

IF(IP.EQ.MP)  GO  TO  22 

21  R(K2)=U2+U3 
R(K5)=U4+U5 

22  K2=K2+N2 
UA=UB 

27  CONTINUE 
12  CONTINUE 

RETURN 
END 
C* ************************** **************SUBROUTINE  ANTFF 
SUBROUTINE  ANTFF (THETA , PHI , THETAS , PHIS , ARAD , ETF , EPF) 
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C  COMPUTES  CLOSED  FORM  ANTENNA  PATTERN  IN  THE  FAR  FIELD. 

C  ETF  AND  EPF  ARE  RETURNED.  ANGLES  ARE  IN  RADIANS. 

C  PRESENT  ENTRY  IS  FOR: 

C 

C  UNIFORM  ILLUMINATION  SCANNED  TO  A  DIRECTION  (THETAS.PHIS) 

C 

COMPLEX  ETF.EPF, JK.SCL.EB 
JK=(0. 0,6. 283185307) 
BK=6. 283185307 
PI=BK/2. 
C  SET  RADIATION  PATTERN  TO  ZERO  IN  THE  REAR  HEMISPHERE 
IF(COS(THETA).LT.O. )  THEN 

EB=(0. ,0. ) 
ELSE 

BB=SIN(THETA)*SIN(PHI)-SIN(THETAS)*SIN(PHIS) 
AA=SIN(THETA)*COS(PHI)-SIN(THETAS)*COS(PHIS) 
ARG=BK*ARAD*SqRT(AA**2+BB**2) 
SCL=-(0. ,1.)*2.*PI**2 
IF(ABS(ARG).LT.1.E-5)THEN 
EB=SCL/2. 
ELSE 

EB=SCL*(BESSJ1(ARG)/ARG) 
ENDIF 
ENDIF 
C  EB  IS  THE  X  COMPONENT;  GET  ETHETA  (ETF)  AND  EPHI  (EPF) 
C  ASSUME  AN  ELEMENT  FACTOR,  ELMT 
ELMT=SQRT(ABS(COS(THETA) ) ) 
ETF=EB*COS (THETA) *COS (PHI ) *ELMT 
EPF=-EB*SIN(PHI ) *ELMT 
RETURN 
END 
i>******#*****«********#****#*****#*****«******BESSEL  FUNCTIONS 

FUNCTION  BESSJ(N.X) 
C  RETURNS  THE  BESSEL  FUNCTION  B  OF  ORDER  N  (>1)  AND  REAL 
C  ARGUMENT  X. 

PARAMETER  (IACC=40,BIGN0=1 .E10 ,BIGNI=1 . E-10) 
IF(N.EQ.O)  THEN 
C  IF  N=0  CALL  BESSJO 

BESSJ=BESSJ0(X) 
ELSE  IF(N.EQ.l)  THEN 
C  IF  N=l  CALL  BESSJ1 

BESSJ=BESSJ1(X) 
ELSE 
C  IF  N>1  USE  RECURSION 
BESSJ=0. 

IF(X.NE.O.)  THEN 
T0X=2./X 

IF(X.GT.FLOAT(N))  THEN 
BJM=BESSJ0(X) 
BJ=BESSJ1(X) 
DO  11  J=1,N-1 
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BJP=J*TOX*BJ-BJM 

BJM=BJ 

BJ=BJP 

11  CONTINUE 
BESSJ=BJ 
ELSE 

M=2*((N+INT(SQRT(FL0AT(IACC*N))))/2) 
BESSJ=0. 

JSUM=0. 

SUM=0. 

BJP=0. 

BJ=1. 

DO  12  J=M,1,-1 

BJM=J*TOX*BJ-BJP 

BJP=BJ 

BJ=BJM 

IF(ABS(BJ) .GT.BIGNO)  THEN 

BJ=BJ*BIGNI 

BJP=BJP*BIGNI 

BESSJ=BESSJ*BIGNI 

SUM=SUM*BIGNI 

ENDIF 

IF(JSUM.NE.O)  SUM=SUM+BJ 

JSUM=1-JSUM 

IF(J.EQ.N)  BESSJ=BJP 

12  CONTINUE 
SUM=2.*SUM-BJ 
BESSJ=BESSJ/SUM 
ENDIF 

ENDIF 

ENDIF 

RETURN 

END 

FUNCTION  BESSJO(X) 
C 

C  BESSEL  FUNCTION  OF  0  ORDER,  REAL  ARGUMENT  X 
C  (SEE  'NUMERICAL  RECIPES',  P. 172) 
C 

REAL*8  Y,Pl,P2,P3,P4,P5,qi,Q2,Q3,Q4,Q5,Rl,R2,R3,R4,R5,R6, 

*  S1,S2,S3,S4,S5,S6 

DATA  PI , P2 , P3  , P4 , P5/ 1 . DO , - . 109862827D-2 , . 27345 10407D-4 , 

*  - . 2073370639D-5 , . 2093887211D-6/ 

DATA  qi ,Q2,Q3,Q4,Q5/-. 1562499995D-1 , . 1430488765D-3 , 

*  -.6911147651D-5,.7621095161D-6,-.934945152D-7/ 

DATA  Rl , R2 , R3 , R4 , R5 , R6/57568490574 . DO , - 1 3362590354 . DO , 

*  651619640 . 7D0 , -11214424 . 18D0 , 77392 . 33017D0 , -184 . 9052456D0/ 
DATA  S1.S2, S3, S4.S5.S6/57568490411. DO, 1029532985. DO, 

*  9494680 . 718D0 , 59272 . 64853D0 , 267 . 8532712D0 , 1 . DO/ 
IFCX.EQ.O.)  THEN 

BESSJ0=1. 
ELSE  IF(ABS(X).LT.8.)  THEN 
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Y=X**2 
BESSJ0=(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6) ) ) ) )/ 

*  (S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6))))) 
ELSE 

AX=ABS(X) 

Z=8./AX 

Y=Z**2 

XX=AX-. 785398164 

BESSJ0=SQRT(.636619772/AX)*(C0S(XX)*(P1+Y*(P2+Y*(P3+ 

*  Y*(P4+Y*P5))))-Z*SIN(XX)*(Ql+Y*(Q2+Y*(q3+ 

*  Y*(Q4+Y*Q5))))) 
ENDIF 

RETURN 

END 

FUNCTION  BESSJl(X) 
C 

C  BESSEL  FUNCTION  B  OF  ORDER  1,  REAL  ARGUMENT  X 
C  (SEE  'NUMERICAL  RECIPES',  P. 173) 
C 

REAL*8  Y,Pl,P2,P3,P4IP5,qi,q2,q3,Q4,qS,Rl,R2,R3,R4,R5,R6) 

*  S1,S2,S3,S4,S5,S6 

DATA  P1,P2,P3,P4,P5/1.D0, . 183105D-2 ,- .3516396496D-4, 

*  .2457520174D-5.-.20337019D-6/ 

DATA  q 1 , q2 , q3 , q4 , q5/ . 04687499995D0 , - . 2002690873D-3 , 

*  . 8449199096D-5 , - . 99228987D-6 , . 105787412D-6/ 

DATA  Rl , R2 , R3 , R4 , RS , R6/723626 14232 . DO , -7895059235 . DO , 

*  242396853 . 1D0 , -2972611 . 439D0 , 15704 . 4826D0 , -30 . 16036606D0/ 
DATA  S 1 , S2 , S3 , S4 , S5 , S6/144725228442 . DO , 2300535 178 . DO , 

*  18583304 . 74D0 , 99447 . 43394D0 , 376 . 9991397D0 , 1 . DO/ 
IFCX.EQ.O.)  THEN 

BESSJ1=0. 
ELSE  IF(ABS(X).LT.8.)  THEN 
Y=X**2 
BESSJ1=X*(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6)))))/ 

*  (S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6))))) 
ELSE 

AX=ABS(X) 

Z=8./AX 

Y=Z**2 

XX=AX-2. 356194491 

BESSJl=SqRT(.636619772/AX)*(C0S(XX)*(Pl+Y*(P2+Y*(P3+ 

*  Y*(P4+Y*P5))))-Z*SIN(XX)*(qi+Y*(q2+Y*(q3+ 

*  Y*(q4+Y*Q5)))))*SIGH(1.,X) 
ENDIF 

RETURN 
END 


C  +  ************  ************************************* SUBROUTINE  CIRCRTP 

C  SUBROUTINE    :        CIRCRTP 

C  DATE         :         1  OCTOBER  1991 
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C  REVISED 
C  PROGRAMMER 
C  2ND  REVISION 
C  PROGRAMMER 


c 

+ant 

+an 

c 

_ 

_ 

c 

i 

\ 

\ 

c 

E      = 

-jn/(4*pi*k)      /_ 

/_ 

c 

x=0 

y=0 

c 

26  February  1992 
R.M.  FRANCIS 

27  October  1992 
K.  A.  KLOPP 

C 

C  COMMENTS  :   THIS  SUBROUTINE  WILL  RIGOROUSLY  CALCULATE  THE  ELECTRIC  FIELD 

C  FOR  A  CIRCULAR  APERTURE.  THE  FIELD  IS  CALCULATED  AT  THE  COORDINATES 

C  SPECIFIED  BY  RH(.)  AND  ZH( . )   THE  APERTURE  IS  LOCATED  AT  Z  =  0 , 

C  AND  SCANNED  TO  A  DIRECTION  (THS.PHS).  THE  SUBROUTINE  OGIVE  IS  THE 

C  SOURCE  OF  THE  GEOMETRIC  DATA  REQUIRED  BY  CIRC  TO  PERFORM  COMPUTATIONS. 

C  ALL  PHYSICAL  DIMENSIONS  ARE  NORMALIZED  TO  WAVELENGTH. 

C 

* 

* 

xl  * 

e  ( )*TAPER(rp)  * 

* 
* 

SUBROUTINE  CIRCRTP (CNPHI , XP , AP , ARAD , THS , PHS , 
*  PHI,RZ,Z,CIRCR,CIRCT,CIRCP) 

INTEGER  CNPHI 

DIMENSION  XP(100),AP(100) 

REAL  K 

COMPLEX  J , JK , Gl , G2 , X 1 , CON , CC , SUMX , SUMY , SUMZ 

COMPLEX  CIRCR.CIRCT.CIRCP 

PI=3. 141592654 

ETA=120.*PI 
C  SET  FIELD  COMPONENTS  TO  ZERO  IN  THE  REAR  HEMISPHERE 

CIRCR=(0. ,0.) 

CIRCT=(0. ,0.) 

CIRCP=(0. ,0.) 
C  Z  OF  THE  OBSERVATION  POINT  IS  >=  0  IN  THE  FORWARD  HEMISPHERE 

IF(Z.GE.O.)  THEN 

SUHX=(0. 0,0.0) 

SUMY=(0. 0,0.0) 

SUMZ=(0. 0,0.0) 

K=6. 283185307 

J  =  (0.0,1.0) 

JK=  (0.0,6.283185307) 

C0N=-J/(4.*PI) 
C  OMIT  THE  FACTOR  OF  ETA  SINCE  IT  IS  NOT  INCLUDED  IN  SUBROUTINE  ZMAT. 
C  ARAD**2  IS  FOR  SCALING  OF  GAUSSIAN  INTEGRATION.   READJUST  SCALING 
C  TO  AGREE  WITH  OLD  SUBROUTINE  SPHERE  (ARBITRARY) 

CC=  CON 
C  OUTER  LOOP:  INTEGRATE  OVER  X,  -ARAD  <  X  <  ARAD. 
C  INNER  LOOP:  INTEGRATE  OVER  Y,  -ARAD  <  Y  <  ARAD. 
C  EVALUATE  ANTENNA  FIELD  AT  RZ.Z: 

XO=RZ*COS(PHI) 

YO=RZ*SIN(PHI) 
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STS=SIN(-THS) 
S=RZ/SQRT(RZ**2+Z**2) 
C=Z/SqRT(RZ**2+Z**2) 
C  INTEGRATE  IN  X,Y  COORDINATES 
DO  50  M=1,CNPHI 
X=ARAD*XP(M) 
DO  60  N=1,CHPHI 
Y=ARAD*XP(N) 
RP=SQRT(X**2+Y**2) 
C  SKIP  INTEGRATION  POINTS  BEYOND  THE  ANTENNA  RADIUS 
IF(RP.LE.ARAD)  THEN 
PHIPRIHE=ATAN2(Y,X+1.E-10) 

R=SQRT(((RZ-RP)**2+Z**2)+(4*RZ*RP*SIN((PHI-PHIPRIME)/2)**2)) 
G1=(((K*R)**2)-1-(JK*R))/R**3 
G2=(3+(3*JK*R)-(K*R)**2)/R**5 
X1=JK*(RP*C0S(PHIPRIME-PHS)*STS-R) 
Y1=X0-X 
X2=Y1**2 
Y2=Y0-Y 

SUMX=SUMX+(CC*AP(M)*AP(N)*(G1+(X2*G2))*CEXP(X1))*TAPER(RP) 
SUMY=SUMY+(CC*AP(M)*AP(N)*Y1*Y2*G2*CEXP(X1))*TAPER(RP) 
SUMZ=SUMZ+(CC*AP(M)*AP(N)*Y1*Z*G2*CEXP(X1))*TAPER(RP) 
ENDIF 
60    CONTINUE 
50    CONTINUE 
C  ASSUME  AN  ELEMENT  FACTOR,  ELMT 
ELMT=SQRT(ABS(C)) 

CIRCT=(C*COS(PHI)*SUMX+C*SIN(PHI)*SUMY-S*SUMZ)*ELMT 
CIRCR=(S*COS(PHI)*SUMX+S*SIN(PHI)*SUMY+C*SUMZ)*ELMT 
CIRCP=(-SIN(PHI)*SUMX+COS(PHI)*SUMY)*ELMT 
ENDIF 
RETURN 
END 
C*******************^**********************FUNCTION  TAPER 

FUNCTION  TAPER(RHO) 
C  SPECIFY  ANTENNA  ILLUMINATION  FUNCTION.   REAL  FUNCTION  OF 
C  RHO  ONLY 

TAPER=1. 

RETURN 

END 
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